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Abstract 

\ 

JA  relative  orbital  element  estimator  is  designed  using 
least  squares  estimation.  Range  and  range  rate  measurements 
are  taken  and  a  vector  of  Delaunay  elements,  relative  to  an 
interceptor  satellite,  is  estimated  using  the  nonlinear 
least  squares  equation.  The  estimator  incorporates  a  state 
vector  coordinate  transformation  from  relative  position  and 
velocity  to  relative  orbital  elements.  Two  modes  of  the 
estimator  program  are  designed,  batch  and  sequential,  where 
the  sequential  mode  uses  Bayes  estimation.  Three  test  cases 
are  analyzed  and  indicate  satisfactory  performance.  Problem 
areas  include  estimator  dependence  upon  a  highly  accurate 
initial  estimate  of  the  element  vector  to  start  the  estima¬ 
tion  process.  The  orbits  are  restricted  to  noncircular  for 
both  satellites  and  the  orbits  must  be  non  coplanar.  Range 
and  range  rate  data  for  the  estimator  are  provided  using  a 
truth  model. 
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DESIGN  OF  AN  ORBITAL  ELEMENT 
ESTIMATOR  USING  RELATIVE  MOTION  DATA 


I .  Introduction 

Tracking  an  orbiting  satellite  from  ground  based  sta¬ 
tions  is  the  primary  method  for  determining  a  satellite's 
set  of  orbital  elements.  In  a  rendezvous  and  intercept 
situation  the  target  orbital  elements  relative  to  the  pri¬ 
mary  or  intercept  satellite  are  needed  to  compute  the 
maneuvers  necessary  to  accomplish  the  intercept.  This  data 
can  be  provided  by  ground-based  tracking  stations  or  can  be 
determined  using  an  estimation  process  on  board  the  inter¬ 
ceptor  that  uses  relative  motion  data  of  the  target  from  the 
interceptor. 

Most  relative  motion  studies  involve  linearization  of 
the  equations  of  motion  about  a  circular  orbit.  Developing 
an  orbital  element  estimator  using  relative  motion  data  is 
based  upon  the  dynamics  or  equations  of  motion  relating  the 
relative  motion  of  two  nearby  satellites.  The  most  commonly 
used  equations  for  relative  motion  are  the  Hill's  equations 
(Ref  5:111).  These  equations  are  used  to  develop  relative 
motion  relationships  between  two  spacecraft  that  are  in 
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neighboring  near-circular  orbits.  This  restriction  to  near 
circular  orbits  limits  the  use  of  Hill's  equations  for  use 
in  an  orbital  element  estimation  process  using  relative 
motion.  Eades  and  Drewry  (Ref  4)  develop  a  set  of  relative 
motion  equations  but  restrict  the  orbits  to  the  same  orbital 
period  and  restrict  one  satellite  to  a  circular  orbit.  The 
estimation  process  used  by  NASA  in  the  Apollo  Guidance 
Computet  Rendezvous  Filter  is  designed  to  estimate  both  the 
position  and  velocity  state  vector  of  the  primary  and  target 
satellite  (Refs  9,  3).  Measurements  used  are  range,  range 
rate,  azimuth,  and  elevation  angles.  There  is  no  restric¬ 
tion  to  the  type  of  orbits,  but  the  Rendezvous  filter  does 
not  estimate  a  set  of  orbital  elements.  A  satellite-to- 
satellite  orbit  determination  system  designed  by  NASA  is  the 
Tracking  and  Data  Relay  Satellite  System  ( TDRSS )  (Ref  6). 
This  system  uses  least  squares  estimation  processing  radar 
measurements  from  a  geostationary,  circular  orbit,  station. 
This  system  is  linked  to  a  ground  station  that  processes  the 
tracking  data  and  therefore  is  not  autonomous.  Also,  the 
tracking  satellite  is  restricted  to  a  circular  orbit. 

Past  studies  of  relative  motion  have  dealt  with  the 
near  circular  orbit  case.  No  method  of  orbital  element 
determination  from  an  in-orbit  platform  using  relative 
motion  for  noncircular  orbit  scenarios  has  been  developed. 
Development  of  a  system  to  determine  orbital  elements  from  a 
spacecraft  has  its  advantages  in  that  it  means  longer  visa- 


bility  and  observation  time  from  space  than  from  the  ground 
tracking  stations.  Determining  the  orbital  elements  is 
necessary  to  plan  and  accomplish  a  rendezvous  and  intercept 
of  two  spacecraft  independent  of  ground  tracking  facilities 
Ground  tracking  data  could  be  supplemented  with  the  data 
determined  by  the  orbiting  tracking  systems  for  selected 
targets. 

The  orbital  element  estimator  designed  in  this  study 
uses  least  squares  estimation  theory  to  estimate  a  set  of 
relative  orbital  elements  of  a  target  from  a  primary 
spacecraft.  Incorporated  in  the  design  is  a  state  vector 
coordinate  transformation  from  relative  motion  data  to  rela 
tive  orbital  elements.  Using  this  transformation  in  the 
least  squares  estimation  process  is  a  unique  concept  in  the 
design  of  this  orbital  element  estimator. 


» 
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II.  Least  Squares  Estimation  Theory 

Least  squares  estimation  is  based  upon  probability 
theory,  specifically  the  principle  of  maximum  likelihood, 
where  the  best  estimate  of  a  state  is  the  value  at  which  the 
probability  of  it  being  the  true  state  is  maximized. 
Probability  Theory 

Equation  (1)  introduces  the  Gaussian  probability  den¬ 
sity  function. 

P  ( e )  =  ((2*)"1/2/o>  exp  (-(e)2/2a2)  (1) 

The  probability  density  function  gives  the  probability 
that  an  error,  denoted  e,  will  lie  between  e  and  6e  as 
P(e)6e.  Integrated  over  the  interval  of  all  possible  error 
this  function  is  normalized  to  unity,  as  shown  in 
Equation  ( 2 ) . 

OO  00 

/  P(e  )de  =  (2ti)~1/2  a"1  /  exp  (-e2/2a?')de  =  1  (2) 

—  OO  —00 

A  gaussian  error  function  is  shown  graphically  in 
Figure  1  and  supplemented  by  Equation  (3). 

P(X)  =  ( 2n ) 1/2  o_1  exp(-(X-XQ)2/2o2)  (3) 


P 
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The  value  XQ  is  the  true  measurement  value.  The 
standard  deviation,  o,  is  a  function  of  the  accuracy  of  the 
measurement  device.  The  Gaussian  probability  of  an  estimate 
being  within  la  of  the  true  value  is  .68  or  68%  probability. 

The  concept  that  links  probability  theory  and  least 
squares  estimation  is  the  principle  of  maximum  likelihood. 
Simply  stated,  the  best  estimate  of  the  true  value  is  the 
estimated  value  which  maximizes  the  probability  of  obtaining 
the  true  measurement. 

The  best  estimate  of  X  is  Substituting  this  for  the 

true  state  yields  the  following  probability  function  for  N 

independent  measurements.  Equation  (4),  where  X  is  the 

S\ 

measurement  vector  and  X  is  the  estimated  measurement 
vector. 

/  N  N  _  «  9  „ 

P(Xi)  =  ( 2n  )  -N/2  p  0  .-1  exp(_r  (Xi-Xi)V2a/)  (4) 

l  1=0  1 

To  maximize  the  probability  the  absolute  value  of  the 
argument  of  the  exponential  is  minimized  using  Equation  (5). 


_d 

A 

dX 


N  _  *  n 

2a .  2 


N  (X.-X^ 

"  =  - 2 

1=1  a . 


=  0 


(5) 


.  .  _  2  . 

Minimizing  (X-X)  is  necessary  to  maximize  the  probabi¬ 
lity,  hence  the  name  least  squares,  where  the  sum  of  the 

_  / \ 

squares  of  the  errors  is  minimized.  (X-X)  is  called  the 

residual,  it  is  the  difference  between  the  observed  measure¬ 
rs 

ment,  X,  and  the  estimate  of  X,  again  denoted  X. 
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Dynamics 

Estimation  of  the  state  of  a  dynamical  system,  such  as 

a  set  of  satellite  orbital  elements,  involves  linearized 

.  /\ 

dynamics.  The  estimate  of  the  state,  X,  is  propagated  in 
time  by  the  vector  differential  shown  in  Equation  (6). 

A  A 

X  =  g(X,  t)  (6) 

For  trajectories  of  the  state  near  the  true  trajectory, 
X=X0+6X,  the  dynamics  of  the  variations  in  the  trajectory  is 
given  in  Equation  (7), 

£  =  ~kQ  +  6*  =  g(XQ  +  6X,t)  (7) 

Using  a  Taylor  series  expansion  yields  Equation  (8). 

£0  +  61  =  g(XQ,t)  +  ?g(X0,t)6X  +  0(2)  (8) 

Equation  (8)  reduces  to  Equation  (9)  and  introduces  the 
A  matrix  which  propagates  the  dynamics  of  the  variations. 


6X  =  A(t) 5X, 


where 

Aij(t)  =  Vg(XQ,t)  = 

The  state  transition  matrix,  $,  propagates  the 
variations  in  the  state  as  shown  in  Equation  (10). 


(9) 


5X(t)  =  $(t,tQ)  6X(t0) 


(10) 
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The  dynamics  of  the  state  transition  matrix  are  deter¬ 
mined  by  Equation  (11). 

*(t,tQ)  =  A(t)*(t,t0)  (11) 

The  covariance,  denoted  P(t),  is  a  matrix  containing  the 
average  squared  errors  of  the  states  related  to  themselves 
and  other  elements.  The  state  transition  matrix  propagates 
the  covariance  by  Equation  (12). 

P(t)  =  *(t,t0)  P  (tQ)  $T  (t,tQ)  (12) 

Linear  estimation  can  be  formulated  using  the 
linearized  dynamics  and  subsequently  related  to  nonlinear 
estimation,  which  is  commonly  used  in  astrodynamical 
problems. 

Linear  Least  Squares  Estimation 

For  a  linear  dynamical  system  the  state  transition 
matrix  propagates  the  state  from  an  epoch  time,  tQ,  to  a 
future  time,  t.  The  state  at  epoch  is  to  be  estimated  from 
data  measurements,  Z^,  taken  over  some  time  interval.  The 
state  vector  is  related  to  the  measurement  vector  by 
Equation  (13)  where  is  the  observation  relation  matrix 
and  e^  represents  the  error  in  the  measurements. 


Associated  with  the  measurement  vector  Z  is  the  Q 


matrix  which  contains  the  squares  of  the  a  values  of  the 
measurement  device. 

Equation  (13)  is  rearranged  to  find  an  expression  for 
the  error  between  measurement  and  the  estimated  measurement, 
the  result  is  Equation  (14). 

ei  =  Zi  -  HiX(t0)  (14) 

The  state  vector  X(t)  is  related  to  the  state  at  epoch 
time  using  the  state  transition  matrix  propagation  X(t0)  to 
X(t).  Incorporating  this  into  Equation  (14)  yields 
Equation  (15). 


ei  =  zA  -#iX(t0)  (15) 

where 

Hi  =  H*<ti,t0) 

Substituting  Qi  as  and  the  error  equation  (Equation  (15), 
into  the  Gaussian  error  function.  Equation  (4),  yields  an 
expression  for  the  probability  of  error  in  the  measurement 
set  (Equation  16). 

P(e)  =  (2tt)  _N/2  |  Q  |  "1/2e(“1/2J)  (16) 

— T  i  — 

where  J  =  e  CT^e. 


To  maximize  the  probability,  J  is  minimized  in  the  same 


way  the  argument  of  the  exponential  was  minimized  in 
Equation  (5). 

A  A 

=  2-s  (z  -  Xo)TQ_1  (Z  -  Xo)  =0  (17) 

ax  ax 

where  X  is  the  estimate  of  the  state  vector 
Differentiating  Equation  (17)  yields  Equation  (18). 


T  -1  -  T  -1- 

ft  Q  ft  x  -  ft-  Q  z  =  0  (18) 

Rearranging  Equation  (18)  leads  to  an  expression  which 
is  the  fundamental  linear  least  squares  equation. 


x(to)  =  (#Tcrtyr1  tq“1z 


(19) 


where  (ft  Q  ft.)  is  the  COvariance  matrix. 


Nonlinear  Least  Squares  Estimation 

Nonlinear  estimation  is  based  upon  linear  estimation 
theory  using  the  linearized  dynamics  and  modified  components 
of  linear  estimation. 

The  measurements  are  now  a  nonlinear  combination  of  the 
states  and  are  given  by  the  observation  relation, 


Z(tt)  =  G(X(ti),(ti))  (20) 

This  relation  is  linearized  to  give  Equation  (21). 
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Hi  is  evaluated  on  the  estimated  or  reference 
trajectory. 

The  residuals  are  given  by  Equation  (22)  using  the 
nonlinear  observation  relation. 

Fi  =  Zi  -  G(Xref ,ti)  (22) 

Following  the  development  of  the  linear  least  squares 

/ 

equation,  using  the  nonlinear  terms  developed,  gives  the 
fundamental  nonlinear  least  squares  equation,  Equation  (23). 

6X(t0)  =  Q_1F  (23) 

and  the  estimated  state  vector  is  given  as  Equation  (24). 

X(tQ)  =  Xref  (tQ)  +  6X(t0)  (24) 

Nonlinear  least  squares  Equations  (23)  and  (24)  are  used 
to  update  the  reference  state  vector  iteratively  until  it 
converges  upon  a  solution.  The  convergence  criteria  is 
based  upon  the  true  solution  being  the  result. 

Theoretically  6X(tQ)  will  converge  to  zero,  in  practice  it 
should  be  allowed  to  converge  to  well  within  the  associated 
square  root  of  the  covariance  for  that  element,  .  The 

residuals,  independently,  should  be  of  order  /OiT*  as  they 
converge. 

The  algorithm  below  shows  the  step  by  step  iterative 
process  used  to  converge  upon  a  solution. 


NONLINEAR  LEAST  SQUARES  ALGORITHM 


1.  From  each  measurement  calculate: 

a.  rt  =  ZL  -  G(Xref  (tiJ'tj.) 

b.  Hi 

c»  $  { t i  / 1^ ) 

2.  Assemble  vector/matrices  necessary  for  nonlinear 
least  squares  equation. 


r  = 

~rl 

r2 

• 

n 

H1  *1*" 

H2  *2 

• 

iO 

II 

( - 

• 

CM 

o 

1— 1 
o 

1 _ 

• 

• 

_ri_ 

• 

• 

-Hi  n- 

1 

• 

• 

o 

H- 

1 

3.  Compute  update  using  least  squares  equation 

6^<  tQ )  =  (^TQ_V)_1  TQ_1r 

4.  Update  reference  solution. 

XUQ)  =  X(tQ)  +  s9(t0) 

5.  Convergence  Check!  If  convergence  criteria  is 

/\ 

met  X{t0)  is  the  solution.  If  not  met,  return  to 
Step  1  where,  using  the  observation  relation,  a 
new  set  of  residuals  are  computed. 

Nonlinear  least  squares  estimation  is  implemented  with 
a  batch  of  measurements  using  the  given  algorithm. 

Combining  sequential  measurements  with  previous  estimates  of 
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a  state  is  also  accomplished  using  the  fundamental  equation. 

A  sequential  nonlinear  least  squares  process  is  imple¬ 
mented  using  the  previous  estimate  of  the  state,  denoted  X- , 
as  data  with  an  associated  covariance,  denoted  P~(t0). 

Since  the  estimate  is  at  epoch  time  the  observation  relation 
matrix  for  this  data  is  the  identity  matrix.  The  components 
of  the  fundamental  equation  for  the  sequential  least  squares 
are  given  in  Equation  (25). 


Substituting  these  into  the  fundamental  equation  for 
nonlinear  estimation  gives  Equation  (26). 

6X  =  (P~(t0)  +^TQ~V)  (p‘(t0)(X  -X)  +^TQ_1?)  (26) 

A 

where  X  is  the  previous  state  estimate. 

Sequential  least  squares  is  more  commonly  called  Bayes 
estimation.  Bayes  estimation  minimizes  the  squares  of  the 
sequential  residuals  and  the  discrepancy  between  the  pre¬ 
vious  estimate  and  the  estimate  of  the  state  vector  deter¬ 
mined  by  the  Bayes  estimation  process.  The  Bayes  estimation 
equation  is  used  iteratively  which  is  similar  to  the  given 
least  squares  algorithm. 


III.  Design  and  Development 

Classical  and  Delaunay  Orbital  Elements 

The  orientation  and  size  of  a  satellite  orbit  is 
described  by  a  set  of  orbital  elements.  The  classical  orbi¬ 
tal  elements  are  the  most  commonly  used.  This  estimator 
uses  the  Delaunay  orbital  elements  which  can  be  directly 
related  to  the  classical  elements.  The  Delaunay  elements 
are  the  standard  set  of  canonical  elements  for  the  two  body 
problem.  All  element  sets  are  mathematically  equivalent, 
and  therefore  one  set  can  be  related  to  another  set. 

To  best  describe  the  Delaunay  elements  the  meaning  of 
the  classical  elements  is  necessary.  Figure  2  gives  the 
classical  orientation  elements. 


Fig  2.  Orientation  Orbital  Elements. 


The  longitude  of  ascending  node  is  SI,  i  is  the  inclination 
to  the  equator  plane,  and  w  is  the  argument  of  perigee.  The 
vernal  equinox  is  the  direction  in  space  used  as  a  reference 
for  the  inertial  coordinate  system  used  in  this  design.  The 
dimensional  elements  specify  the  size  and  shape  of  the  orbit 
and  relate  position  in  the  orbit  with  time.  The  semi-major 
axis  is  denoted  a,  e  is  the  eccentricity,  it  defines  the 
shape  of  the  orbit  and  M  is  the  mean  anomaly  which  deter¬ 
mines  the  satellite’s  position  in  the  orbit.  Mean  anomaly 
is  an  angular  element.  Commonly  used  with  the  classical 
elements  is  a  coordinate  system  that  has  unit  vector  ‘p’ 
directed  to  perigee,  the  closest  point  in  the  orbit  to  the 
attracting  body,  'q'  is  perpendicular  to  ’p'  and  completes 
the  orthogonal  right  handed  set.  This  frame  is  called  the 
perifocal  coordinate  system,  or  PQW,  as  shown  in  Figure  3. 


Apogee 


W 


Fig  3.  Perifocal  Coordinate  System. 
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The  Delaunay  elements  can  be  related  to  the  classical 
elements  using  Equations  (1)  through  (6)  (Ref  11:24). 

L  =  / pa,  where  p  is  the  gravitational  constant  (1) 


1  “  n(t-TPERIGEE)  =  M'  Where  "  =V5"  (2> 
G  =  /y a(l-e2)  =  L  J 1-e2  (3) 
g  =  w  (4) 
H  =  /pa( 1-e cos  i  =  G  cos  i  (5) 
h  =  ft  (6) 


Delaunay  elements  g  and  h  are  identical  to  the  classi¬ 
cal  orientation  elements  w  and  ft.  Element  1  is  the  mean 
anomaly  and  is  also  defined  by  Equation  (7)  which  introdu¬ 
ces  E,  the  eccentric  anomaly,  an  angular  measurement  from 
the  center  of  the  ellipse. 

1  =  E  -  e  sin  E  (7) 

Element  L  is  related  to  the  semi-major  axis  and  used  to 
define  the  total  energy  of  the  orbit.  G  is  related  to  the 
eccentricity  and  is  the  total  angular  momentum  of  the  orbit. 
H  is  the  inertial  vertical  component  of  the  angular 


momentum 


The  classical  orbital  elements  in  terms  of  the  Delaunay 


elements  are  given  in  Equations  (8)  through  (13). 


a  =  L2/u  (8) 
e  =  /l  -  G2/L2  (9) 
i  =  cos-1 ( H/G ) ,  sin-1 ( 1-H2/G2 )  (10) 
a)  =  g  (11) 
U  =  h  (12) 
M  =  1  (13) 


Equations  of  Motion  and  Relative  Motion 

Relative  motion  can  be  described  as  the  apparent  motion 
of  another  object  as  seen  from  the  primary  object.  The 
relative  motion  of  two  nearby  satellites  is  found  by  deter¬ 
mining  the  position  and  velocity  of  each  vehicle  and  dif¬ 
ferencing  the  vectors  accordingly  to  find  the  relative  posi¬ 
tion  and  velocity.  The  position  and  velocity  must  be  in  a 
common  coordinate  system.  The  position  and  velocity  of  a 
satellite  is  a  function  of  the  orbital  elements  and  time. 
Figures  4  and  5  show  how  the  relative  position  and  velocity 
vectors  are  found.  (Figures  4  and  5  follow  on  the  next 
page.  ) 

Position  and  velocity  computations  for  each  satellite 
are  initially  computed  in  the  PQW  coordinate  system.  To 
find  the  inertial  relative  motion  the  vectors  are  trans¬ 
formed  into  an  inertial  frame,  defined  IJK,  where  I  is 
directed  to  the  vernal  equinox.  Figure  6  shows  the  rela¬ 
tionship  between  PQW  and  IJK  coordinate  frames. 
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Interceptor 


Lg  5  . 


Relative  Velocity  Commutation. 


The  mean  motion  is  defined  as  dl/dt  and  is  given  in 
Equation  (4). 


dl 

dt 


2  /T  3 
V  /L 


(4) 


dE/dl  is  the  change  in  the  eccentric  anomaly  with 
respect  to  1,  the  mean  anomaly.  Differentiating  implicitly 
yields  Equation  (5). 

dl  =  dE  -  e  cos ( E )dE  (5) 


Substituting  e  in  terms  of  the  Delaunay  elements  in 
Equation  (6) , 

dl  =  dE  -  /l-G2/L2  cos  E  dE  (6) 

and  separating  leads  to  Eq  (7) 

dE  _  1 

dT - — 7 - 2 — 7 — - 

1  -  /l-GVL  cos  E 

dr/dE  is  found  by  differentiating  the  vector  components 
of  position  with  respect  to  the  eccentric  anomaly  E. 
Equations  (8)  and  (9)  are  the  results. 


(8) 


drp/dE  =  (-L2/y)  sin(E) 
drg/dE  =  (LG/ u)  cos(E) 


(9) 


f 


Substituting  expressions  from  Equations  (4),  (7),  (8), 
and  (9)  into  (3)  gives  expressions  for  the  velocity  vector 
in  the  PQW  frame.  Equations  (10)  and  (11). 

VP  =  -u  sin( E)/( L( 1  -  /l-G2/L2  cos ( E ) ) )  (10) 

Vq  =  pG  cos  (  E )/(  L2  ( 1  -  A-G2,/L2  cos  ( E ) )  )  (11) 


The  position  and  velocity  state  vector  is  given  as  a 
six-component  vector  Equation  (12). 


X  PQW  _ 


rp“ 

rQ 

rw 

Vp 

VQ 


,  where  rw  and  Vw  =  0 


LVWJ 


Rotation  from  PQW  to  IJK  frame  is  accomplished  using 
Equation  (13),  where  R,  the  transformation  matrix  from  PQW 
to  IJK  is  given  in  Appendix  A. 


XIJK  =  [r]Xpqw  (13) 

The  position  and  velocity  of  each  satellite  is  computed 
at  given  times  and  differentiated  accordingly  in  the  iner¬ 
tial  frame  to  yield  the  relative  position  and  velocity 
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vector  of  the  target  from  the  interceptor  (Equation  14). 


XRELATIVE 


rT , I  "  rI,l’ 

■*rl' 

rT,J  -  *1'J 

*rJ 

rT,K  “  rI,K 

= 

ArK 

VT,I  -  VIfI 

AVj 

VT,J  ~  VIfJ 

AVj 

VT,K  "  VI,K 

avk 
»  • 

(14) 


T  =  Target 
I  =  Interceptor 

Range  and  Range  Rate  Measurements 

Range  and  range  rate  measurements  are  used  as  obser¬ 
vational  data  for  the  estimator.  These  measurements  were 
selected  since  they  relate  directly  to  relative  position  and 
velocity  and  are  typically  measured  by  a  radar  unit.  The 
radar  unit  used  in  the  design  has  a  100  meter  standard 
deviation,  a,  in  range  accuracy  and  a  0.3  meter  per  second 
o  in  range  rate  accuracy.  A  radar's  accuracy  is  dependent 
upon  several  factors.  Range  to  the  target  and  power  pro¬ 
vided  to  the  radar  are  basically  those  factors  that  affect 
the  accuracy  and  capability.  Cross  section  of  the  target 
and  shape  are  also  important  factors  to  consider.  This 
estimator  uses  a  theoretical  radar  that  is  not  range 
restricted  and  has  the  constant  a  values  given  previously. 

It  is  assumed  that  the  interceptor  radar  is  able  to  track 
any  target,  regardless  of  size,  shape,  or  range.  Special 
cases  tested  used  a  range  restriction  of  less  than  200  Km, 
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which  is  a  reasonable  rang.e  limit  for  a  radar.  Generally, 
the  range  was  not  restricted  for  testing  the  performance  and 
accuracy  of  the  estimator. 

Range  is  measured  along  a  line  of  sight  vector  from  the 
interceptor  to  the  target.  Range  is  the  magnitude  of  this 
line  of  sight  vector,  as  shown  in  Figure  7. 

Target 


Fig  7.  Range  Computation. 

In  terms  of  the  relative  position  the  range  is  given  in 
Equation  (1). 

Range  =  (arj2  +  Arj2  ArK2)^2  (1) 


Range  rate  is  the  rate  of  change  of  the  range  scalar. 
The  dot  product  of  the  rela.tive  velocity  vector  and  the 
relative  position  vector  divided  by  range  gives  range  rate, 
as  shown  in  Figure  8. 


Range  Rate  =  (V  •  r)/|r|  =  | V j  cose 

Fig  8.  Range  Rate  Computation. 

Equation  (2)  gives  the  range  rate  in  terms  of  the  rela¬ 
tive  position  and  velocity. 

Range  rate  =  (AVjArj  +  AVjArj  +  AVKArK)/Range  (2) 

The  estimator  program  reads  the  range  in  distance  units 
(DU),  where  1  DU  =  6378.145  Km,  range  rate  in  distance  units 
per  time  unit  (DU/TU),  where  1  DU/TU  =  7.905376  Km/sec,  and 
time  in  TU,  where  1  TU  =  806.8136  seconds.  The  time  used  is 


the  time  past  epoch,  the  epoch  time  being  when  the  intercep¬ 
tor  was  a  perigee.  (Ref.  1:429) 

Truth  Model 

The  truth  model  generates  the  true  range  and  range  rate 
data  at  given  times  for  various  orbital  scenarios.  The 
truth  model  data  is  generated  using  a  computer  program  that 
propagates  the  position  and  velocity  of  two  satellites,  an 
interceptor  and  target,  in  two  orbits.  The  position  and 
velocity  are  propagated  using  the  equations  previously 
developed,  Equations  (1),  (2),  (10),  and  (11)  of  the 
Equations  of  Motion  section  of  this  chapter.  Time  is  the 
common  variable  that  is  synchronized  for  both  satellites. 
Each  satellite  is  started  from  an  initial  position  and  moved 
using  1  minute  time  increments.  At  each  time  the  associated 
eccentric  anomaly  is  computed  for  the  equation  of  motion. 
Orbit  size  and  orientation  are  computed  using  a  given  set  of 
orbital  elements  for  each  satellite.  The  inertial  position 
and  velocity  vector  at  each  measurement  time  are  computed 
and  differenced  accordingly  to  give  the  relative  position 
and  velocity.  Range  and  range  rate  are  computed  using  the 
relative  position  and  velocity  and  are  given  by 
Equations  (1)  and  (2),  which  are  repeated  from . the  Range  and 
Range  Rate  section. 

Range  =  (Arj2  +  Arj2  +  Ar  k2)"^2  (1) 

Range  Rate  =  (AVjArj  +  AVjArj  +  AVRArK)/Range  (2) 
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The  truth  model  computer  program  gives  range,  range 
rate,  and  time  at  10  minute  intervals  for  an  orbit  scenario, 
but  minor  modifications  can  adjust  the  time  interval.  The 
computer  program  and  associated  flow  chart  for  the  truth 
model  are  presented  in  Appendix  C. 

Relative  Orbital  Element  Estimator 

The  estimator  designed  in  this  study  uses  least  squares 
estimation  theory  to  estimate  a  set  of  relative  orbital  ele¬ 
ments.  The  relative  orbital  element  state,  Equation  1,  is 
the  difference  between  the  known  interceptor  elements  and 
the  estimated  target  elements  at  an  epoch  time,  here  chosen 
to  be  the  time  when  the  interceptor  is  at  perigee  prior  to 
taking  range  and  range  rate  measurements.  Initially  an 
estimate  of  the  relative  element  vector  is  needed  to  start 
the  estimation  process.  In  this  study  the  initial  estimate 
is  found  by  perturbing  the  true  relative  element  vector  by  a 
certain  percentage  selected  by  the  user. 
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The  estimator  is  designed  using  the  developed  theory, 
however,  a  modification  unique  to  the  design  is  incor¬ 
porated.  The  modification  used  in  this  estimator  that  is 
different  from  typical  least  squares  estimation  is  the  state 
vector  coordinate  transformation  from  relative  position  and 
velocity  to  a  relative  orbital  element  state  vector. 
Variations  in  the  relative  position  and  velocity  are  related 
to  the  variations  in  the  Delaunay  orbital  elements.  This 
modification  is  incorporated  via  the  state  transition 
matrix. 

From  the  Least  Squares  Theory  section,  Equation  (2) 
depicts  the  classical  way  the  4>  matrix  is  used  to  propagate 
the  variations  of  the  state  vector  in  time. 

5X(t)  =  ♦(t,t0)5X(t0)  (2) 

The  $  used  in  this  design  is  not  classical.  Here 
4>  propagates  the  variations  of  the  position  and  velocity 
relative  to  the  interceptor  orbit.  These  variations  are 
propagated  from  an  epoch  time  to  a  particular  time,  in  this 
case  a  measurement  time.  As  well  as  propagating  motion, 

*  also  relates  the  variations  in  position  and  velocity  to 
variations  in  the  Delaunay  elements.  Thus,  *,  called  the 
state  transition/state  vector  coordinate  transformation 
matrix,  relates  relative  position  and  velocity  to  the  rela¬ 
tive  orbital  elements.  $  is  formed  by  developing  rela¬ 
tionships  between  the  satellite  position  and  velocity  in  the 
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PQW  coordinate  frame  and  the  associated  Delaunay  orbital 
elements  and  then  transforming  the  result  into  the  IJK  coor¬ 
dinate  frame.  The  position  and  velocity  vector  (Equation  3) 
and  the  orbital  elements  vector  (Equation  4)  are  related  to 
each  other  in  the  PQW  reference  frame  by  the  relation  given 
in  Equation  (5) 

Xr#v  =  [rp,  Tq,  rw,  Vp,  Vq,  Vyf]  (3) 

^ELEMENTS  =  hJ  (4) 


Equation  (5)  results  in  a  6x6  matrix  relating 
variations  of  the  six  components  of  position  and  velocity  to 
variations  of  the  six  orbital  elements  in  the  PQW  frame. 

The  result  of  Equation  (5)  is  transformed  into  the  IJK  coor¬ 
dinate  frame  using  Equation  (6). 


d(Xr,V) IJK 

d<W 


+  [R] 


3*xelm^  r,v 


d<Xr,V> 

d(XELEM) 


R  is  a  coordinate  rotation  matrix  from  the  PQW  to  IJK 
frame.  The  result  of  Equation  (6)  is  a  6x6  matrix  relating 
the  variations  of  relative  position  and  velocity  to 
variations  in  the  Delaunay  elements  in  the  IJK  frame.  The 
detailed  vector  and  matrix  components  for  the  development  of 
$  as  well  as  the  other  components  necessary  for  implemen- 
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tation  of  the  nonlinear  least  squares  equation  are  pre¬ 
sented  in  Appendix  A. 

The  validity  and  accuracy  of  4  is  checked  using  numeri¬ 
cal  derivatives  and  comparing  them  with  the  appropriate 
column  vector  of  <J>  at  the  given  time.  A  position  and  velo¬ 
city  vector  is  determined  using  a  given  set  of  orbital  ele¬ 
ments.  One  element  is  perturbed  a  small  amount,  5,  and 
another  position  and  velocity  vector  is  found.  The  vector 
difference  divided  by  the  change  of  the  element  should  agree 
with  the  associated  column  vector  of  4.  Equation  (7) 
depicts  the  check  method. 

— — —  -  i§  (7, 

ELEMENT  ELEMENT 

After  4  was  formulated  several  checks  using  various  orbit 
types  indicated  the  validity  and  accuracy  of  $. 

Incorporating  the  state  vector  coordinate  transfor¬ 
mation  into  this  estimator  meant  that  relative  position  and 
velocity  data,  easily  measured,  could  be  directly  trans¬ 
formed  to  the  orbital  elements  necessary  for  orbital  deter¬ 
mination.  The  estimator  computer  program  and  flowchart  are 
presented  in  Appendix  D. 

The  estimator  is  checked  using  data  generated  by  the 
truth  model.  Various  scenarios  are  simulated  where  range 
and  range  rate  measurements  of  a  target  satellite  are  made 
from  an  interceptor  satellite  at  selected  time  intervals.  A 
batch  mode  and  a  sequential  mode  are  available  where  an 
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estimation  of  the  relative  orbital  elements  is  made  based  on 


a  single  batch  of  measurements  or  several  sequential  batches 
taken  in  time  segments.  The  sequential  mode  uses  Bayes 
estimation  theory.  Flexibility  in  the  truth  model,  its  abi¬ 
lity  to  generate  data  for  any  interceptor-target  scenario, 
leads  to  flexibility  in  the  ways  the  estimator  can  be 
tested,  therefore  revealing  strong  points  or  weaknesses  in 
the  accuracy  and  versatility  of  the  design.  The  Testing  and 
Results  section  of  this  report  presents  the  statistics  and 
discusses  the  various  scenarios  and  special  tests  performed 
using  the  relative  orbital  element  estimator  program. 
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IV.  Testing  and  Results 


The  performance  of  the  estimator  in  the  test  cases  is 
evaluated  by  comparing  the  statistics  of  the  estimated  ele¬ 
ment  vector  to  the  statistics  predicted  by  the  estimator. 

The  covariance  matrix  for  the  batch  least 

squares  and  (P{-)-^  for  the  Bayes  or  Sequential 

estimation  give  the  estimators  prediction  of  how  well  it  can 
estimate  the  state  vector  or  relative  orbital  elements.  The 
diagonal  elements  of  the  covariance  matrix  give  the  varian¬ 
ces,  <j2,  for  each  element.  These  variances  are  indicators 
of  the  estimator's  predicted  performance.  From  the  esti¬ 
mated  relative  orbital  elements  resulting  from  a  set  of 
simulations,  where  random  noise  corrupted  the  measurements , 
a  sample  variance  is  computed  using  Equation  (1),  which  is 
the  definition  of  a  sample  variance. 

2  l  N  _  _  2 

°e  =  N^T  J=1  (XJ,e  ESTIMATED  ”  Xe ,TRUE  ^  (1) 

where  _ 

X  is  the  relative  element  state  vector 
N  is  the  number  of  simulations 
e  denotes  the  Delaunay  Element 

The  variances  computed  using  Equation  (1)  are  the 
results  of  the  estimations  made  in  the  simulations  for  a 
particular  scenario.  Similarity  between  these  actual  ele¬ 
ment  variances  and  the  variances  predicted  by  the  estimator 
from  the  covariance  matrix  indicate  satisfactory  performance 
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of  the  estimator. 


The  measurement  device  accuracy  over  all  ranges  to  the 
target  is  given  in  Equations  (2)  and  (3). 

GRANGE  =  57X10_5  yy  =  200  meters  (2) 

0 RANGE  RATE  =  3.85X10“^  DU/TU  =  0.3  meters/second  (3) 

These  values  indicate  the  accuracy  limits  to  which  the 
measurement  device  operates.  From  the  residuals  an  indica¬ 
tion  of  the  accuracy  to  which  the  estimator  updated  the  pre¬ 
dicted  measurement  can  be  found  and  compared  to  the  given 
radar  accuracies. 

The  root  mean  squared  value  (RMS)  of  the  residuals  are 
computed  in  each  simulation  using  Equation  (4),  the  defini¬ 
tion  for  RMS  value  for  a  sample  group  of  measurements. 

Measurement  RMS  =  i  (rm  -  rp)2  (4) 

where 

rm  is  the  measurement 

rp  is  the  estimator  predicted  measurement 

N  is  the  number  of  measurements  taken. 

The  variances  for  the  range  and  range  rate  residuals 
are  computed  using  Equation  (5)  using  the  measurement  RMS 
values  from  Equation  4. 
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where 


N-r  -  nrNo 

Nr  is  the  number  of  simulations 

No  is  the  number  of  observations  per  simulation 

The  square  root  of  the  variance  gives  the  standard 
deviation.  The  standard  deviation  of  the  residuals,  com¬ 
puted  using  Equation  (5),  are  compared  to  the  given  standard 
deviation  of  the  radar  accuracies.  Agreement  of  the  asso¬ 
ciated  standard  deviations  indicates  the  estimator  has  pre¬ 
dicted  the  range  and  range  rate  measurements  to  the  accuracy 
of  the  measurement  device. 

Three  test  case  scenarios  are  statistically  evaluated 
to  check  the  performance  of  the  estimator.  Simulations  of 
each  test  case  are  made  using  random  noise  inputs  to  the 
measurement  data.  The  initial  estimate  used  to  start  the 
estimation  process  is  the  true  relative  element  vector  for 
that  particular  test.  The  true  vector  is  used  because  dif¬ 
ficulties  are  experienced  when  a  perturbed  initial  element 
vector  is  used  to  start  the  estimator.  These  difficulties 
in  the  requirement  of  a  highly  accurate  initial  estimate  are 
discussed  shortly  in  this  section.  The  test  cases  are  now 
described . 

Case  I  uses  the  least  squares  batch  mode  of  the 
estimator.  In  this  mode  a  number  of  range  and  range  rate 
measurements  are  taken  and  processed  to  give  estimate  of  the 
relative  element  vector.  Case  I  uses  13  observations  over  a 
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1  orbit  period  of  the  interceptor.  Table  1  gives  Case  I 
classical  orbital  elements  and  relative  Delaunay  elements. 


TABLE  1 

CASE  I  ORBITAL  ELEMENTS 

Interceptor  Target 


a 

1.3  DU 

1.35 

e 

0.2 

0.3 

i 

45° 

46° 

n 

45° 

46° 

w 

45° 

46° 

lo 

0° 

0° 

A  L 

0.2171958 

DU2/TU 

Al 

0  RAD 

AG 

-0.00876201 

du2/tu 

Ag 

.0174533 

RAD 

AH 

-.0199932 

du2/tu 

Ah 

. 0174533 

RAD 

1Q  is  the  mean  anomaly  at  the  epoch  time.  The  epoch 
time  is  selected  as  the  time  at  which  the  relative  orbital 
elements  are  to  be  estimated.  Figure  9  shows  the  orbit  con¬ 
figuration  at  epoch  time  for  Case  I.  Case  II  also  uses  the 
least  squares  batch  mode.  The  primary  difference  is  that  in 
Case  II  at  epoch  time  there  is  a  0.1  radian  difference  in 
the  mean  anomaly.  In  Case  I  at  epoch  the  two  satellites  are 
co-linear  on  the  P  vector  of  the  PQW  coordinate  frame  while 
in  Case  II  there  is  an  out  of  phase  or  non  co-linear 
situation  at  epoch  as  shown  in  Figure  10.  Case  II  elements 
are  identical  to  Case  I  except  the  target  1Q  is  0.1  radian 
so  the  Al  at  each  epoch  is  0.1  radian. 
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Case  III  uses  the  sequential  mode  of  the  estimator.  In 
this  mode  a  least  squares  batch  estimation  based  on  an  ini¬ 
tial  set  of  observations  is  made.  This  is  an  initial  set  of 
13  observations.  This  is  followed  by  sequential  sets  of 
observations  where  two  sets  of  three  measurements  over  a  30 
minute  time  segment  are  made.  After  the  initial  set  and  two 
sequential  sets  of  measurements  are  taken,  the  relative 
orbital  elements  are  estimated.  Case  III  elements  are  iden¬ 
tical  to  Case  II. 

Appendix  B  presents  the  measurement  data  for  each  test 
case.  This  data  is  obtained  from  the  truth  model,  in  the 
simulations  random  noise  corrupts  the  measurements.  An 
explanation  of  the  random  noise  inputs  to  the  data  is  also 
contained  in  Appendix  B. 

Statistics  for  Case  I  are  based  on  25  simulations,  Case 
II  is  based  on  30,  and  Case  III  is  calculated  from  18  simu¬ 
lations.  In  each  simulation  different  sequences  of  random 
noise  inputs  are  used.  As  noted  earlier  in  this  section, 
the  estimator  uses  the  true  relative  element  vector  in  all 
evaluated  test  cases  as  an  initial  estimate  to  start  the 
estimation  process,  therefore,  the  test  cases  presented 
indicate  the  estimator  error  in  estimating  a  solution  for 
the  relative  element  state  vector,  when  given  the  true  state 
vector.  The  statistics  for  the  test  cases  are  included  in 
Tables  2  through  7,  shown  on  the  following  pages.  Contained 
in  the  first  column  of  Tables  2,  3,  and  4  are  the  variances 
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computed  from  Equation  (1)  using  the  resultant  vectors  from 
the  associated  case  simulations.  The  second  column  contains 
the  average  variances  for  each  element  from  the  covariance 
matrix  for  that  case.  Similarly  in  the  first  column  of 
Tables  5,  6,  and  7  is  the  standard  deviation  computed  using 
Equation  (5)  and  the  second  column  contains  the  given  para¬ 
meters  for  the  radar. 

The  comparison  of  statistics  indicates  good  agreement 
between  the  actual  and  predicted  statistics.  The  estimator 
predicted  statistics  agree  with  the  statistics  of  the  resi¬ 
duals  and  estimated  element  vector.  This  agreement  tells  us 
that  the  estimator  is  functioning  optimally  for  the  given 

scenarios  and  initial  conditions.  Comparing  magnitudes  of 

2  2 

the  variances  show  that  a^l.  and  °GG  are  smaller  than  the 
other  element  variances.  This  indicates  the  estimator's 
ability  to  determine  the  dimensional  element  L  and  G  with 
more  confidence  or  accuracy  than  the  angular  elements  1,  g, 
h,  H.  In  comparing  Case  II  and  III  variances,  where  the 
conditions  at  epoch  are  the  same,  but  Case  II  uses  the  batch 
mode  while  Case  III  incorporates  more  data  using  the  sequen¬ 
tial  mode,  the  variances  for  Case  III  are  notably  smaller 
than  Case  II.  This  decrease  in  the  covariance  from  the 
batch  to  sequential  case  indicates  the  advantage  of  incor¬ 
porating  more  data  sequentially.  A  decrease  in  the 
covariance  implies  more  confidence  in  the  estimator's  pre- 
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TABLE  4 


CASE  III  VARIANCE  COMPARISON 


Element 

From  Element  Vector 
°ii2 

Equation  (1) 

From  Covariance  Matrix 

Oii2 

L 

1.19  E-12 

3.2  E-13 

1 

1.09  E-10 

1.5  E-10 

G 

4.5  E-12 

2.3  E-12 

g 

1.44  E-8 

4.05  E-9 

H 

9.58  E-9 

1.6  E-9 

h 

2.69  E-9 

7.7  E-9 

TABLE  5 

CASE  I  MEASUREMENT  STANDARD  DEVIATION  COMPARISON 


Measurement 


°R,RR  .  °R, RR 

Computed  from  Residuals  Radar  Parameters 

_ Equation  (5) _ _ 


Range  1.96  E-5  DU  1.57  E-5 

Range  Rate  3.96  E-5  DU/TU  3.85  E-5 


TABLE  6 


CASE  II  MEASUREMENT  STANDARD  DEVIATION  COMPARISON 


Measurements 

°R,  RR 

Computed  from  Residuals 
Equation  (5) 

°R,  RR 

Radar  Parameters 

Range 

1.74  E-5 

1.57  E-5 

Range  Rate 

3.77  E-5 

3.85  E-5 

TABLE  7 

CASE  III  MEASUREMENT  STANDARD  DEVIATION 

COMPARISON 

Measurement 

°R,RR 

Computed  from  Residuals 
Equation  (5) 

°R,  RR 

Radar  Parameters 

Range 

1.55  E-5 

1.57  E-5 

Range  Rate 

3.53  E-5 

3.85  E-5 
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diction  of  the  element  set  when  the  Bayes  estimator  process 
is  implemented. 

Simulations  using  initial  estimates  other  than  the  true 
vector  were  also  tested.  These  initial  estimates  were  com¬ 
puted  using  a  perturbing  factor  that  perturbed  the  true  vec¬ 
tor  by  a  percentage  determined  by  the  magnitude  of  the  per¬ 
turbing  factor.  The  sign  of  the  perturbation,  +  or  -,  is 
determined  by  a  random  process.  These  simulations  were  made 
using  noise  free  measurements  and  the  initial  estimate  of 
the  relative  element  vector  was  perturbed.  The  percentage 
of  perturbation  to  the  true  element  vector  was  limited  to 
values  less  than  or  equal  to  1%.  When  values  greater  than 
1%  were  used  the  estimator  would  not  converge  to  a  solution. 
When  noise  corrupted  measurements  were  used  the  limit  of  the 
perturbing  percentage  decreased  to  approximately  1/10%. 

This  indicated  that  with  noise  corrupted  measurements  the 
estimator  is  unable  to  converge  to  a  solution  unless  a 
highly  accurate  initial  estimate  of  the  element  vector  is 
provided . 

Several  orbit  scenario  cases  indicated  weaknesses  in 
the  estimator  performance.  If  the  target  and  interceptor 
orbits  are  coplanar  the  residuals  will  always  diverge  and 
convergence  to  a  solution  is  impossible.  The  relative 
inclination  in  this  case  is  zero,  therefore  the  relative 
longitude  of  ascending  node,  h,  is  undefined.  The  equation 
for  satellite  dynamics  and  reference  frame  transformations 
contain  h  and  in  the  coplanar  case  with  h  undefined  a  singu- 
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larity  exists.  This  singularity  occurs  for  low  relative 
inclination,  near  zero,  situations.  For  circular  orbits 
perigee  is  undefined  and  in  near  circular  orbits  the  perigee 
may  be  difficult  to  locate  due  to  small  changes  in  the 
distance  from  the  attracting  body  throughout  the  orbit 
trajectory.  Since  the  estimator  uses  a  perigee  dependent 
reference  coordinate  system  the  circular  orbit  weakness  was 
recognized.  Since  the  classical  and  Delaunay  elements  are 
dependent  upon  the  existence  of  a  perigee  the  estimator  has 
difficulty  estimating  the  angular  elements  1  and  g  that  need 
a  defined  perigee  for  definition.  Simulations  revealed  that 
the  estimator  could  not  estimate  1  and  g  separately.  This 
was  indicated  by  large  variances  for  the  angular  elements; 
however,  the  sum  of  1  and  g  gave  an  accurate  estimation  of 
target  angular  position  at  epoch  time.  Therefore,  in  near 
circular  cases  determining  angular  position  in  the  orbit  is 
possible  but  determining  the  angular  elements  1  and  g 
separately  results  in  erroneous  values  for  those  individual 
elements.  The  estimator  estimates  Delaunay  elements  L  and  G 
more  accurately  than  the  angular  elements  1,  g,  h,  and  H  as 
indicated  by  the  estimator  predicted  variances  found  in 
Tables  2  through  7. 
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f.  Conclusions  and  Recommendations 

Least  squares  estimation  theory  and  techniques  were 
used  in  the  design  of  the  relative  orbital  element 
estimator.  The  new  and  unique  feature  incorporated  in  this 
estimator  is  the  state  vector  coordinate  transformation  from 
relative  position  and  velocity  coordinates  to  relative  orbi¬ 
tal  elements.  It  is  possible  that  the  difficulties  encoun¬ 
tered  in  the  orbit  scenarios  discussed  in  the  results  can 
be  solved  using  several  modifications  or  changes  to  the 
development  and  implementation  of  the  design.  Several  solu¬ 
tions  or  recommendations  are  now  presented. 

Using  the  conventional  set  of  Delaunay  elements  lead  to 
some  of  the  difficulties  experienced  with  circular,  near 
circular,  and  coplanar  orbits  of  the  interceptor  and  target. 
In  a  circular  orbit  perigee  is  undefined.  With  perigee 
undefined,  the  argument  of  perigee,  g,  is  therefore  unde¬ 
fined.  Since  the  equations  of  motion  and  coordinate  trans¬ 
formation  from  the  perifocal  to  inertial  coordinate  system 
are  dependent  upon  g  as  a  variable  these  equations  become 
incalculable.  In  near  circular  orbits  the  element  g  is 
difficult  to  define  since  the  orbit  is  so  near  being 
circular  and  the  problems  of  singularities  in  the  circular 
orbit  case  arise.  Coplanar  orbits  have  the  relative  incli¬ 
nation  equal  to  zero.  In  this  situation,  longitude  of 
ascending  node,  h,  is  undefined.  Again  equations  for  motion 
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and  transformation  will  contain  an  undefined  term  and  be 
incalculable. 

These  problems  can  be  solved  using  another  set  of  orbi¬ 
tal  elements  that  can  successfully  handle  low  inclination, 
circular  or  near  circular  orbits.  The  ideal  set  for  this  is 
the  Equinoctal  orbital  elements.  These  elements  are  used  by 
the  NORAD  Space  Computational  Center  (Ref  2)  and  are  free 
from  both  zero  inclination  and  eccentricity  difficulties. 

The  Equinoctal  elements  are  given  in  Equations  (1)  through 


(6). 

(1) 

a  £  =  e  cos  it 

where  u=g+w  is  the  longitude 

(2) 

ag  =  e  sin  n 

of  perigee 

(3) 

M  =  mean  motion 

(4) 

L=ft  +  <i)  +  M  =  ir 

+  M  =  mean  longitude  when  M  is  the 

mean  anomaly 

(5) 

sini  sing 
x  1  +  cosi 

(6) 

sini  cosn 
^  1  +  cosi 

With 

circular  orbit  elements  af  and  ag  are  zero,  no  elements 

are 

undefined  in  this 

case.  For  zero  inclination  elements 

X  and  are  zero  and  again  no  elements  are  undefined. 
Reformulating  the  state  transition/coordinate  state  trans¬ 
formation  using  the  equinoctal  elements  is  recommended  to 
eliminate  the  coplanar  and  circular/near  circular  orbit 


difficulties. 


The  equation  for  position  and  velocity  for  the 
satellite  is  initially  computed  using  a  perifocal  coordinate 
system.  These  calculations  are  dependent  upon  an  orbit 
where  a  perigee  exists.  Selection  of  a  coordinate  system 
not  dependent  upon  the  existence  of  a  perigee  will  eliminate 
this  problem.  All  computations  in  an  inertial  coordinate 
system  is  one  way  of  solving  this  problem. 

The  difficulties  experienced  when  the  initial  estimate 
of  the  relative  element  vector  is  perturbed  from  the  true 
relative  element  vector  are  of  primary  concern  since  these 
difficulties  would  limit  the  scope  with  which  the  estimator 
could  be  used  operationally.  A  good  assumption  is  that  the 
initial  estimate  used  to  start  the  operation  of  the  estima¬ 
tor  will  come  from  a  ground-based  tracking  facility  and  not 
be  the  true  set  of  orbital  elements  due  to  limitations  in 
the  tracking  accuracy.  Starting  with  the  initial  estimate 
the  interceptor,  using  range  and  range  rate  measurements, 
can  estimate  a  more  accurate  set  of  orbital  elements  for  the 
target.  Presently  a  highly  accurate  initial  estimate  is 
required  to  insure  convergence  upon  a  solution.  This  highly 
accurate  initial  estimate  may  in  fact  be  more  accurate  than 
the  set  the  estimator  could  provide.  Therefore,  the  use  of 
the  estimator  would  introduce  more  error  to  the  orbit  deter¬ 
mination  process.  Reformulating  or  modifying  the  state 
transition/coordinate  transformation  matrix  using  the 
Equinoctal  elements  may  increase  estimator  versatility. 
Adding  relative  angular  measurements  to  -supplement  the  range 


and  range  rate  measurements  may  also  enhance  the  estimator's 
ability  to  successfully  converge  to  a  solution  starting  with 
a  less  accurate  initial  estimate  of  the  element  vector. 
Incorporating  more  types  of  observational  data,  such  as 
angular  measurements,  would  give  more  information  to  the 
estimator  to  process.  Angular  measurement  may  also  increase 
the  accuracy  of  the  estimator's  prediction  of  angular 
elements,  thereby  decreasing  the  variances  of  the  angular 
elements  to  values  comparable  to  the  dimensional  elements. 

Modifying  relative  position  and  velocity  estimators, 
such  as  the  NASA  Apollo  Rendezvous  Filter,  with  the  state 
vector  coordinate  transformation  technique  developed  in  this 
design  can  create  a  new  orbit  determination  method  to  be 
used  by  satellite  tracking  facilities  from  a  space-based 
platform,  such  as  the  Space  Shuttle  or  an  orbiting  station. 
The  modified  estimator  could  be  designed  and  tested  simi¬ 
larly  to  the  method  used  in  this  study.  The  truth  model 
would  have  to  be  modified  to  give  relative  angular  measure¬ 
ments  and  a  coordinate  transformation  from  inertial  to 
vehicle  reference  frame  would  be  required. 

Due  to  the  requirement  of  a  highly  accurate  initial 
estimate  necessary  to  start  the  estimation  process  this 
design  is  not  yet  considered  by  me  to  be  operational. 
However,  with  continued  research  and  development  of  the 
ideas  and  techniques  used  in  this  study  along  with  implemen¬ 
tation  of  my  recommendations  can  bring  us  closer  to  develop- 
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ment  of  an  operational  space  based  orbital  estimation  pro¬ 
cess  to  supplement  and  enhance  our  ground  based  tracking 
facilities  in  tracking  and  cataloging  the  many  orbiting 
vehicles  in  space,  both  friendly  and  others. 


V 
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Appendix  A 

Vector/Matrix  Components  of  the  Least  Squares 
Estimation  Equation 

1.  State  Transition/State  Coordinate  Transformation 
Matrix,  $ 

The  $  used  in  this  design  is  not  a  classical  4>  matrix. 
*  relates  the  relative  position  and  velocity  state  vector 
and  the  relative  orbital  element  state  vector.  $  also  pro¬ 
pagates  the  state  vector  in  time.  $  is  formed  by  develop¬ 
ing  relationships  between  the  satellite  position  and 
velocity  in  the  perifocal,  PQW,  coordinate  frame  and  the 
Delaunay  orbital  elements  and  then  transforming  the  result 


into  the 

inertial,  IJK, 

coordinate 

frame . 

Position  and 

velocity 

are  denoted 

here  as  Xrv 

and  the 

orbital 

elements 

are  denoted  X^lm  an<3 

are 

defined 

in 

Equations  (1) 

and  (2) 

Xrv  = 

lrI' 

rJ'  rK ' 

Vi 

'  Vj, 

vK] 

(1) 

XELM  = 

[L, 

lr  Gf  gf 

H, 

h] 

(2) 

Equation  (3)  gives  the  PQW  position  and  velocity  rela¬ 
tion  to  the  elements. 


d ( *rv } PQW 
d(^ELM)  6x6 
Transforming 


d ( Xr, v ] PQW  dE  +  3<Xr,v}PQW  (3) 

dE  ’  d<W  3<W 

Equation  (3)  into  the  IJK  frame  is  accomplished 


by  Equation  (4). 


—  — 

-  - 

d<*r,v>IJK 

_  P  3R  1 

(X  )  +  [R] 

...  r ,  v  PQW 

d(Xr,v)PQW 

d(XELM) 

L3<wJ 

(4) 


where 


R  =  Rotation  matrix  PQW  to  IJK  frame. 

$  is  a  6x6  matrix  given  as  Equation  (5). 


<t>  = 


d(Xr#v,IJK 


d(XELM) 


(5) 


The  rotation  matrix  is  a  function  of  g,  h,  G,  H.  R  is 
given  in  Equation  (6)  where  C  denotes  cos  and  S  denotes  sin 
function.  R  transforms  a  position  or  velocity  vector  from 
the  PQW  to  IJK  coordinate  frame. 


[R]  = 


C(h)C(g)-S(h)S(g)  (H/G)  -C ( h )  S (g ) -S ( h ) C (g  )  ( H/G )  S(h)A-H2/G2  ' 
S(h)C(g)+C(h)S(g) (H/G)  -S ( h ) S ( g ) +C ( h ) C (g ) ( H/G )  -C(h)  A-H2/G2 


S(g) 


/l-H2 


C(g )  A-H2/G2 


H/G 


(Equation  6) 


Equations  (7),  (8),  (9),  and  (10)  give  the  partial  deriva¬ 
tive  of  R  with  respect  to  the  Delaunay  elements. 

f-C(h) S(g ) +S(h )C(g ) ( H/G)  S( h ) S (g ) H/G-C (h )C(g )  ( 


3  R" 


3  R~ 

ag_ 


-S(h)S(g)+C(h)C(g) (H/G)  -S ( h ) C (g ) +C ( h ) S ( g ) ( H/G )  0 

C(g )  /l-H2/G2  -S(g)  A-H2/G2  0 

(Equation  7) 


49 


1 F  ■  1  "JT*  •”  1,11 


f 


-S(h)C(g)+(H/G)S(g)C(h)  S(h)S(g)-C(h)C(g) ( H/G)  C(h)/l-H2/G2 
C(h)C(g)-S(h)S(g) ( H/G )  -C ( h ) S ( g ) +S ( h ) C (g ) ( H/G )  S(h)/l-H2/G2 

0  0  0 
(Equation  8) 


S(h)S(g) (H/G2) 

-C ( h ) S ( g ) H/G2 

S (g ) { H2/G3 )//l-H2/G2 
(Equation  9) 


S ( h ) C ( g ) H/G2  S(h) (H2/G3)//l-H2/G2 

-C { h ) C ( g ) H/G2  -C ( h ) ( H2/G3 ) /  /l-H2/G2 
C (g ) (H2/G3)/l-H2/G2  -H/G2 


-S(h)S(g)/G 

C(h)S(g)/G 


-S (g ) ( H/G2 )//  1-H2/G2 
(Equation  10) 


-S(h)C(g )/G 
C(h)C(g )/G 

-C (g  )  (  H/G2  )/ A-H2/G 


-S(h)  (H/G2)/A-H2/G 
C(h)  (H/G2)/A-H2/G2 


1/G 


2 


The  derivatives  of  the  position  and  velocity  in  the  PQW 
frame  with  respect  to  the  eccentric  anomaly  are  given  in 
Equations  (11)  through  (16). 

drp/dE  =  ( -L2  jj  )  sin  (  E-E0 )  (11) 

drQ/dE  =  (LG/m)  cos(E-Eq)  (12) 

drw/dE  =0  (13) 


pL[sin2(E)/ 1-G2/L2  -cos ( E) ( l-/l-G2/L2cos ( E) ) ] 
dVp/dE  =  -  (14) 

[L(l-/1-G2/L2  cos ( E ) ) ] 2 


mi  Ai 
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p 


d  Vq/<3  E  = 


yiL2G[cos(E)sin(E)  /l-G2/L2  -  (l-/l-G2/L2cos(  E)  )  sin  (  E )  )  ] 


[L2 ( 1- /l-G2/L2  cos ( E ) ) ] 2 


dV, 


W 


dE 


=  0 


(15) 

(16) 


The  partials  of  the  position  and  velocity  in  the  PQW  frame 
with  respect  to  the  Delaunay  Elements  are  given  in 
Equations  (17)  through  (29). 


3r 


3  L 
3r 


p  _ 


—  cos(E) 


-  it  A-gJl2  -  ^  (7 . \ . 2)(G2/L3) 

w  p  /l-G  /ir 


3L 

3r 


Q  _  G 


—  sin( E) 


W 


3  L 

3VI 
3  L 


=  0 


(17) 

(18) 

(19) 


ySin(E) [l-/l-G2/L~  COS ( E )  -(G2/L2 )cos( E)//T-G2/L2 ] 

I L ( l-/l-G2/L2  co s ( E ) ) ] 2  (20) 


3  V  _  _______ 

Q  vCOS(E) [ L2 { 1-/ 1-G2/L2  cos ( E) ) -cos ( E)G2//1-  G2/L2) 


3G  _ 

(L2  ( 1-/ 1-G2/L2  cos ( E) ) ] 2 

(27) 


3  V 
W 

-  =0  (28) 

3G 


3  X 

PQW 


The  derivative  of  E  with  respect  to  L  is  given  in 
Equation  (30). 


dE  /  2  2  2  3  /  2  2 

-  =  [(1/Zl-G  /L  )(G  /L  )sin(E)]/[l~/l-G  /L  cos  ( E) )  (30) 

dL 


The  derivative  of  E  with  respect  to  1  is  given  in 
Equation  (31). 
dE  ~  ]  2  7 

—  =  U/(l-/l-G  /L  cos  (  E )  ) )  (31) 

dl 


The  derivative  of  E  with  respect  to  G  is  given  in 
Equation  (32). 


dE  /  2  T  2  /  2  2 

—  =  ((1//1-G  /L  )G  sin(E)/L  ]/(l-/l-G  /L  cos(E)]  (32) 

dG 
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2.  Observation  Vector  and  Relation 

Range  and  range  rate  are  measured  quantities,  and  each 
time  these  measurements  are  taken  the  following  vectors  and 
matrices  are  assembled.  The  observation  vector  Z,  is  given 
in  Equation  ( 34 ) . 


Z 


Range  to  Target 
Range  Rate 


(33) 


Z  is  a  function  of  the  relative  position  and  velocity, 
as  shown  in  Equation  (34). 

(Arj.2  +  Arj2  +  ArK2)"^ 

z  =  (AVIArl  +  AVjArj  +  AVR ArK )/RANGE  (34) 

If  the  relative  position  and  velocity  are  defined  in 
vector  XjJK  then  Z  is  given  as  Equation  (35)  where  the 
observation  vector  is  a  nonlinear  function  called  the  obser¬ 
vation  relation  G. 


Z  =  G( X , t ) 


(35) 


Relating  the  state  vector  XIJK  to  the  observation  vec¬ 
tor  Z  is  accomplished  using  Equation  (36). 

H  =  3G(^,t-)-  (36) 

3  X 

In  detail  H  is  given  in  Equation  (37). 
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H 


a RANGE  a RANGE  3 RANGE 

aAr-j.  3Arj  .  3AVK 

a RANGE  RATE  3 RANGE  RATE 

aArj.  3AVk 


(37) 


where  the  components  of  H  are  given  in  Equations  (38) 


through  (40).  Note  that  range  is  not  a  function  of  relative 
velocity. 

3  RANGE  _  ArI,J,K  nft) 

3 Arr  _  „  RANGE 

L  I  d  r  IS. 


a RANGE  RATE 
3Ari,J,K 

3  RANGE  RATE 

a  a  v. 


(RANGE)AVt  ,  „—  (  A  V  •  Ar  )  Ar  _  T  ,,/RANGE 
1  /  J  /  1\  1/  J  f  K 


RANGE' 


Ar 


I,J,K 


RANGE 


I,J,K 

fj-l  is  given  by  Equation  (41). 

y  =  H  $ 

so  in  detail^/  is  represented  by  Eq  (42). 


y- 


d RANGE 


dRANGE 


dL  dl 

dRANGE  RATE 
L  DL 


dRANGE 

iJH 


dRANGE  RATE 
dh  . 


(39) 

(40) 


(41) 


(42) 


3 .  Q  Matrix 

The  Q  matrix  is  a  diagnol  matrix  of  the  values  for 
the  measurement  device.  1  orange  =  100  meters  and 
1  ORANGE  RATE  =  .3  meters/second  were  used  for  the  estimator 
radar.  The  size  of  the  Q  matrix  is  determined  by  the 
number  of  measurements  taken,  it  is  a  2Nx2N  matrix  where  N 
is  the  number  of  measurements  taken  (Equation  43). 
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Residual  Vector 


A  residual  is  defined  as  the  difference  between  the 
measured  value  and  the  estimation  of  the  measurement 
Equation  ( 44 ) . 


rl,N  “ 


where 

X  is  the  relative  position  and  velocity  vector. 

5.  Nonlinear  Least  Squares  Equation 

When  the  components  thus  far  computed  are  substituted 
into  the  nonlinear  least  squares  equation.  Equation  (45), 
the  result  is  Equation  (46),  correction  to  the  reference 


RANGE  1 


RANGE  RATE  1 


RANGE  N 
RANGE  RATE  N 


-  |c(xlfN,tlfN)J 


Test  Case  Data 

Test  case  data  is  generated  using  the  truth  model  com 

puter  program  contained  in  Appendix  C. 

Test  Case  1^  Data 
Classical  Orbital  Elements 


Element  Interceptor  Target 

a  1.3  DU  1.35  DU 

e  0.2  0.3 

i  45°  46° 

w  45°  46° 

Si  45°  46° 

M0  0°  0° 


Delaunay  Orbital  Elements 
Target  -  Interceptor 

Value 

0.02171957876309  DU2/TU 
0°  RAD 

0.008762011386843  DU2/TU 
0.01745329252222  RAD 
0.01999321219953  DU2/TU 
0.01745329252222  RAD 

Tracking  Target  From  Interceptor  Over  1  Orbit  at  10  Minute 
Intervals,  No  Noise  Corruption  On  Data 


Range 

Range  Rate 

Time 

(DU) 

( DU/TU ) 

(TU) 

(MIN) 

0.17513109 

0.09800249 

0.74366638 

10 

0.21953387 

0.02029584 

1.4873328 

20 

0.211625212 

-0.03500555 

2. 23099914 

30 

0.17893223 

-0.04364209 

2.97466552 

40 

0.16317941 

0.01209514 

3.71833189 

50 

0.20293606 

0.09071542 

4.46199827 

60 

0.28797655 

0.13177235 

5.20566465 

70 

0.39208177 

0.14477618 

5.94933103 

80 

0.49886571 

0.13946822 

6.69299741 

90 

0.59442894 

0.11326908 

7.43666379 

100 

0.65887875 

0.05266586 

8.18033017 

110 

0.65981660 

-0.05861977 

8.92399655 

120 

0.567448959 

-0.18229518 

9.66766293 

130 

Conversions 

6378145  m 

=  1  DU 

7.0905376  m/s 

=  1  DU/TU 

806.8136 

=  1  TU 

A  Element 

AL 

Al 

AG 

Ag 

AH 

Ah 
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Test  Case  2  Data 
Classical  Orbital  Elements 


Element  Interceptor  Target 

a  1.3  DU  1.35  DU 

e  0.2  0.3 

i  45°  46° 

u  45°  46° 

a  45°  46° 

Mc  0  0.1  RAD 


Delaunay  Orbital  Elements 
Target  -  Interceptor 

A  Element  Value 

AL  0.02171957876309  DU2/TU 

Al  0.1  RAD 

AG  -0.008762011386843  DU2/TU 

Ag  0.01745329252222  RAD 

AH  -0.01999321219953  DU2/TU 

Ah  0.01745329252222  RAD 

Tracking  Target  From  Interceptor  Over  1  Orbit  at  10  Minute 
Intervals  Between  Data,  No  Noise  On  Data 


Range 

Range  Rate 

Time 

(DU) 

(DU/TU) 

(TU) 

(MIN) 

0.33780912 

0.08567198 

0.74366638 

10 

0.36540682 

-0.00667171 

1.48733276 

20 

0.337044592 

-0.06346114 

2.23099914 

30 

0.27953501 

-0.08556216 

2.97466552 

40 

0.21942767 

-0.06804477 

3.71833189 

50 

0.192346278 

0.00356536 

4.46199827 

60 

0.22684532 

0.083173361 

5.20566465 

70 

0.30400591 

0.117095657 

5.94933103 

80 

0.392584098 

0.116469891 

6.69299741 

90 

0.47041812 

0.08780420 

7.43666379 

100 

0.514020451 

0.02200318 

8.18033017 

110 

0.491439874 

-0.08888904 

8.92399655 

120 

0.384083659 

-0.18700137 

9.66766293 

130 

Conversions 

6378145  m 

=  1  DU 

7.0905376  m/s 

=  1  DU/TU 

806.8136 

=  1  TU 
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Test  Case  2  Data 
Classical  Orbital  Elements 


Element 


Interceptor 

1.3  DU 
0.2 
45° 

45° 

45° 

0 


Target 

1.35  DU 

0.3 

46° 

46° 

46° 

0.1  RAD 


Delaunay  Orbital  Elements 
Target  -  Interceptor 

A  Element 


Value 


AL  0.02171957876309  DU2/TU 

Al  0.1  RAD 

AG  -0.008762011386843  DU2/TU 

Ag  0.01745329252222  RAD 

AH  -0.01999321219953  DU2/TU 

Ah  0.01745329252222  RAD 

Tracking  Target  From  Interceptor  Over  1  Orbit  at  10  Minute 
Intervals  Between  Data.  Two  Batches  of  Three  Sequential 
Measurements  Follow  Initial  Set 


Range 

Range  Rate 

Time 

(DU) 

(DU/TU) 

(TU) 

(MIN) 

Batch  1 

0.24317394 

-0.173231099 

10.4113293 

140 

0.1431328 

-0.08826717 

11.15499568 

150 

0.124736178 

+0.041287686 

11.89866206 

160 

Batch  2 

0.18910942 

0.11915443 

12.64232844 

170 

0.29104299 

0.15135906 

13.38599482 

180 

0.41097969 

0.16973908 

14.1296612 

190 

Noise  Corruption 

Noise  corruption  to  the-  data  is  determined  by  random 
number  generator  that  determines  the  position  along  a  square 
noise  function  modelled  with  a  la  value  equivalent  to  a 
Gaussian  o  value.  The  square  function  is  superimposed  over 
the  Gaussian  distribution.  The  equivalent  lo  value  is  found 
by  equating  the  expressions  for  for  the  gaussian  and 
square  probability  functions  and  finding  the  square  a  in 
terms  of  the  Gaussian  function  a.  This  is  shown  below. 


(1)  a 


square 


-W 


/ 


w 


x2  ^  dx 
2W 


Square 


Gaussian 


Gaussian 


/2'rf  t 


2  2 

e"x  'o! 

dx 


1  3  ,  1, 

3  X  2W 


W 

-W 


2vr 

"5w 


w 

1 


Gaussian 

So  the  square  function  a,  called  W,  is  equivalent  to 
/To  Gaussian,  Equation  (2). 

(2)  W  =  ^Square  =  ^3"°Gaussian 
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Appendix  C 

Truth  Mod e 1  Flowchart  and  Computer  Program 


Enter  Classical  Orbital 
Elements  and  M0,(l0) 

For  Interceptor  and  Target 
Enter  Time  Unit 


Convert  to  Delaunay 
Elements 


* Kepler  Algorithm 


Compute  Position  and  Velocity 
of  Target  in  PQW  Frame,  Then 
Rotate  to  Inertial  Frame 


6j 


TRUTH  MODEL  COMPUTER  PROGRAM 
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Appendix  D 


Relative  Orbital  Element  Estimator  Flowchart 
and  Computer  Program 


I 

i 


Real  Variables 
AX 

ARGP 

CK 

C1-C6 

DXPDE 

DXQDE 

DXWDE 

DVPDE 

DVQDE 

DVWDE 

DEDL 

DEDLS 

DEDG 

DLSCDE 

DL 

DLS 

DG,  DGS ,  DH,  DHS 
E,  ECEN,  ECENT 
ECC 


DEFINITIONS  OF  VARIABLES  USED 
IN  MAIN  PROGRAM 


Serai  major  axis  of  interceptor  (DU) 

Argument  of  perigee  of  interceptor 
(deg,  RAD) 

Difference  between  the  computed  and 
actual  mean  anomaly  in  the  Kepler 
Algorithm 

Diagnol  elements  of  the  Covariance  matrix 

Drp/DE,  derivative  of  the  position  vector 
component  along  the  perigee  direction 
with  respect  to  eccentric  anomaly 

DrQ/DE 

Dr^j/DE 

DVp/DE 

DVq/DE 

DVw/DE 

DE/DL,  Derivative  of  the  eccentric 
anomaly  with  respect  to  the  Delaunay 
element  L 

DE/Dl 

DE/DG 

Dl/DE ,  used  in  Kepler  algorithm 

Ltarget  -  ^interceptor 
■'•target  “  '■interceptor 
See  above 


ECCTGT 


Eccentric  anomaly 
Eccentricity  interceptor 
Eccentricity  target 


G 

GI 

GS 

GSI 

H,  HI,  HS, 

INCL 

L 

LI 

LS,  LST 
LSI 
LSC 
LNOD 

MEANMO 

MU 

NOSR 

NOSRR 

PTB 

PI 

RANGE 

RRATE 

RANG EG 

RRATEG 

RKM 

RRKM 

SIGR 


Delaunay  element  G  target 

Delaunay  element  G  interceptor 

Delaunay  elemeng  g  target 

Delaunay  element  g  interceptor 

HSI  See  above 

Inclination  of  interceptor  orbit 

Element  L  for  target 

Element  L  for  interceptor 

Element  1  target 

Element  1  interceptor 

1  computed  in  Kepler  algorithm 

Longitude  of  ascending  node  of 
interceptor 

Mean  motion  of  orbit 

Gravitational  constant,  1  DU-^/TU2 

Noise  for  range 

Noise  for  range  rate 

Pertubation  input  for  elements 

3.141592654 

Range 

Range  rate 

Range  computed  by  estimator 
Range  rate  computed  by  estimator 
Range  converted  to  KM 
Range  rate  converted  to  KM/sec 
Gaussian  a  for  range  measurement 
Gaussian  a  for  range  rate  measurement 


SIGRR 


TIME 


Time  in  TU 


TIMEM 

Time  in  minutes 

VPPL 

3VP/3L 

VQPL 

3Vq/3L 

VWPL 

3  3  L 

VPPG 

3  Vp/ 3  G 

VQPG 

3Vq/3G 

XPPL 

3  r  p/ 3  L 

XQPL 

3  Tq/ 3  L 

XWPL 

3rw/3L 

XPPG 

3  r  p/ 3  G 

XQPG 

3  Tq/ 3  G 

Integer  Variables 

B1 

Number  of  measurements  x  2 

B2 

i 

' — 1 

CD 

CIO 

Counter 

COUNT 

Counter 

CYCLE  ... 

Counter 

IER 

Error  parameter  for  MATRIX  inversion 
subroutine 

J 

Counter 

K 

Counter 

ID,  ID1 

Measurement  of  identifiers 

MODE 

Mode  of  program,  1  if  batch, 

0  if  sequential 

NUMMEA 

Number  of  measurements  to  be  taken 

NOSE 


Noise  parameter,  1  if  noise  corruption 
to  data 


N 

NN 

SEED 

STORE 

RR 

XX 


Counter 

Counter 

Start  parameter  for  random  number 
generator 

Counter 

Counter 

Counter 


Arrays  (Vectors  and  Matrices ) 


PRPG 

PRPGS 

PRPH 

PRPHS 

XPQW 

R 

DXDL 

DXDLS 

DXDG 

XIDL 

XIDLS 

XIDG 

I  DOS 


[3R/3G]6x6 
[  3  R/  3  g  ]  g  x  g 
[ 3  R/ 3  H] gx6 
[ 3R/3h] gxg 

Position  and  velocity  vector  in  PQW 
frame,  XgxX 

Rotation  matrix  PQW  -*•  IJK,  [ R]  6x6 

(dX/dL]6xl 

[dX/dl]  6xl 

IdX/dG]  6xl 

tdXINERTIAL/dLl 6x1  XINERTI AL  is  x 
rotated  to  IJK  frame 

[dX j/dl ] gxX 
[dXj/dG]  gxX 
[dXx/dg]6xl 
IdXj/dH]  6xl 
(dXj/dhl6xl 


<  I  DM 


TRANS 


T 


$  matrix,  composed  of  column  vectors 
dXj/dL,l,G,g ,H,h.  This  is  the  state 

transition/coordinate  transformation 
matrix  $ 


DEL 

Vector  of  the  differences  between  the 
target  and  interceptor  elements 
[AL,  Al,  AG,  Ag,  AH,  ^^^6x1 

STATE 

Vector  of  relative  position  and  velocity 
[Arj,  Arj,  ArK,  AVj ,  AVj,  AVK]6xl 

HCURL 

^matrix  for  measurement  [^]t2x6 

DELDEL 

Update  vector  for  relative  element  state 
vector,  computed  by  estimator 

WKAREA,  WK 

Work  matrix  for  inversion  subroutine 

ROTT 

Rotation  matrix  PQW  +  IJK,  3x3 

ROT  I 

Rotation  matrix 

POS 

Position  vector  in  PQW  frame 

VEL 

Velocity  vector  in  PQW  frame 

HH 

Matrix  of  for  all  measurements 

l^l'^2'  #3  *  *  •  2nx6 

HHTP 

HH  Transpose 

QINV 

Q  matrix  inverse 

HHTPQ 

Q  used  in  least  squares  equation 

HTPQH 

Res 

Residual  vector 

HTPQH I 

t^Q-1*/)-1 

PHTP 

PHTPQ 

(tfT  q-W-V’Q'I 

POS  I 

Position  vector  in  inertial  frame 

VELI 

Velocity  vector  in  inertial  frame 

POST 

Position  vector  of  target  in  PQW  frame 
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VELT 

Velocity  vector  of  target  in  PQW  frame 

POSTI 

Position  vector  of  target  in  inertial 
frame 

VELTI 

Velocity  vector  of  target  in  inertial 
frame 

HMAT 

H  matrix  for  measurement  time 

RNG 

Vector  of  range  measurements 

RRT 

Vector  of  range  rate  measurements 

POSS 

Storage  matrix  for  position  vectors 

VELS 

Storage  matrix  for  velocity  vectors 

TSTOR 

Vector  of  time  of  measurement  vectors 

PM 

Initial  inverse  covariance  matrix 
computed  for  Bayes  filter,  P— 1 ( — ) 

DELM 

Initial  vector  of  differences  of  elements 
for  Bayes  filter,  5(-) 

HCURLT 

HT 

The  following  matrices/vectors  are  used  in  the  Bayes  filter 
the  main  program 

Ml 

JTq-1 

PP 

P(-)  ♦JTq-1 

PP1 

(P(-)  +^TQ-y)-l 

M3 

f/TQ"1r 

DELDIF 

S(-)  -  ^ESTIMATED 

M4 

P  (  -  )  ~^  (  <5  ( -  )  -  5eST  ) 

TOT 

M3  +  m4 

DELDEL 

PP1  x  TOT,  DELDEL  is  the  Bayes  filter 
computed  updates  to  the  relative  element 
state  vector 

CC 

Vector  composed  of  the  square  roots  of 
the  diagnol  elements  of  the  Bayes 
computed  covariance  matrix 
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RELATIVE  ORBITAL  ELEMENT  ESTIMATOR  PROGRAM 


=  kOG:AI  Si3rIL 


set*  five  it  E1E  N”  PrOG-lkK 


peal  l,ls,g,gs  ,h,  hs,  :,  -i-  ,rrc,  oyfc  t  ,3*30-;  .cy-c^lv/fle,  ;n::, 

«OV'^;  ,  Yf  =  l  ,  <3r-L  ,  X  -fFL  ,7  "r'-  C'***  L,  y-,ct  ,'■£'-•_  ,  3E3 Li  »*  .- PC-,  *3 3 G,  V  =  =  G, 
•t/OPG.OiXG,  <  »  '  »»  L  ‘.Ci?  « ‘  G  =  ,:!  EL,  ’  *7  ‘  ,-..t.G£,.  Ul,  ,  C.GE  G,  £  =  A7EG, 

♦ecctgt,  -  ?c  cr  :t;T,siv  , : ~j c, , o' : , usi ,  , g± , h : » !  e •  c  3 , .  5 r ,  ls c , . n  » 
«DL5C?t,c',  ;i,c'2*: ;  .c  , r-  ,'■*  <i».-.  :s-.,  * 35 =p, 

*gl  ,01:  ,!  ls  ,oi' .,:  G,o‘.i  r-  Er.cGGi ,  j-*,oui,dws,:h;.:  ,  -  t- 
iM-iGei  =: ,v'.i ■  .1: : ,  e2,y v,  ir, t  *cli  ,  j,<,  < , , i _ ^ • ■(. ,c:  jmi  , :i '  ,  stops. , 
*'40S-:,ir:,“.G?,Gf'3 

VJL  =  -  p;-  (  ,,  ; >  ,  =  •  =H  C  '->  ,cr.f  uS(t  ,  r  ) ,=  ’.FGl-  (.'•, 7  )  ,XFCW<  o,  1>  ,P  <--,  i>  , 

•  oxcl  ( 1  ,i> ,  ci  <  ,:>, x:ot  3 <■>,;.>  ,x;  :lc  ,:),  1  c  ep  c  ,  1 ) , 

»xi3Ga,:>  .n^st:  ,:>  ,x::v'<-  ,n  ,x*r  ~3( ,,:), re*. c  <  , 

•3r.Lt.  ,d  ,  s : -  <£ » .  > ,  t:ir  r  ,-. ) ,  1.1:  £■-<.,  1) ,  ■<  <  a  cl  i.ort;,:), 

•POT  I  (  3  .  )r;SCill|J:'LCi;)|Hl<(-  ,0) ,  -IS  =  I  ,  0  i  t<  V  t  .  ■  ,  V  • ) 

»,  mh:  to  (t  , :  ps -tc  (  , :  > ,  ‘Tpoh;  (  if  j)  ,s  -i !  f  t  j,  *  i , 

*=H7P3(t,  .  ) ,  -  C  S*  t :  » 1 ) »  t*7  L  7  ( •: ,  1 )  ,  A  j  ST  t  7 , 1 )  ,  V :  L  i  t ,  * ) ,  =  3  >  7 1  l  •  ,  i  ) , 
•'/SU  £  (  2,1)  ,  P-AT  t  E,  '•)  ,  £  I”-  '*1  ,t  i£.  V  t;  ,o)  ,K  ;G  t£  .  >  ,-  E  I  <2  . 1  , 

♦pcs'  t  ’,2  ) ,  j:.i sc-,  r. ),  ’S"r:  t:  )  ,7s  t  a,  3),i:£Li  it- ,i )  ,u:  jp.t  t o,2> , 

*i2  (•  ,  c )  ,cf  t  j,  •  ) , :  »•..), 

»g£l;:p c  *i  >,t . T(  - , :> ,  •  >r  ,D  ,c:  < . ) 

CO~"-.T.  /.L.1/  l,.  i  »  G  ,3  %',-!.,£,  ECO  , -HJ 
CO'CT  =■ 
tit=: . 


5TC-  -'  = 

CYCLED 
:i  = 

=  I  =  2.1i  l:  S  £  i  ' 

5IG-  7 

SIGFP  =  2.CE  -  ■ 

ap.lS]  .  Y  =  U  I  *  C  CE  D c  IF  lEATSTT  SSCU.C7;,:  .  F  3SY.S  * 

PEAT  * ,*or  . 

st at i  v£ :t;=  :•  is  c?£- : n- ~  is  the  isctg-  cp  aelai.v.  -lexeme 

3„,I‘*r  ;y  =  !:T  ‘CIS-  =  '  ..r  =  i  7>ttN  yOIIc  Is  133£D  T3  "  A  S  U 

* = £w£rTS  • 

^SAC  4  ,*  OS 

opi'.r  •,»  :mP"T  A;1:'-."-  N't*'?  £•  Gt  M  '•  '•£  SEES  ' 

=  £  AC  *  ,  = 

CALL  f  ,'lA‘i  <£'. -C) 

p pit  t  •• ,«  7cs  •••  ■■  tr.i::,L  t  -:agu'<- ieyts  r;  r • 

P£AD  ,!.'!►  -£A 

=  ■■  It  T  •  .J  '  ££=  ;j  •£ 

3  l  =  t  J f<  *  :  A  •  £ 

go  ■  :=i,r. 

OC  A  J=l,31 
=  1 1  V  ( I  i  I )  =  . 

COt.7 1  M’£ 

CO'l'it.J-. 

32=31-: 

301  1=. .  2,~ 

°i v t : , i )  =  t  .  ca»  c 
con*; ‘it.1'. 

30  2  1  =  2,  ■  1 .  u 

3 1  w  1 1 ,  i )  =  .  r  ■  *  3 

CONY! SU£ 

=  (.:••*  -y-  •  si-*  -rj-'  a»is  ci'), .  »c.ty  ,LcyGirJ3£  3-*  • 

• .  •  y  »r  a;  :  'f  £  ;  *-.•  j-.G3£ESi  * 


z; 


ssrt-cppi/n  . 

HS=L  *  OC'PJ/l^  . 
L=SORT<rX) 

g=l*  ((!.-£::•  -  z)  *  ■  .c) 

IdCL  =  If-CL*  Pr/1‘  . 

H  =  G'CC'  <!NCL> 

call  rct eicctij 

5SI=3S 

HSI=HS 


L  I  =L 

3 1  =G 

HI=H 

<X  =  -1 

PH  K!  I 

■ , •  INPUT 

PRINT 

.*  .ST-.-C 

READ  •' 

,CL,Cwr,GG 

pf'  IliT 

- , *  TRUTH 

t  OTt' 

<  ni  i  •  r 

1  L  f  >  ^ 

3  RIM* 

•  ,ch,  •  • , r> 

pp.InT 

* , *  ett ::  ■ 

REAR  • 

>FTF 

3  N I M 

• ,*  -T3  = 

eft  J?t :::  s  T  i  ■-  ■)■•<  - 1:  -  l 


;<sRi  Kr  (  > 

i f  .it.  • . : )  r Mi i 

3 ai:  1 1)  =  r^*  =  : <•• 

ELSE 

j-niiii'a*"?1 

?No:r 

CK=‘  I'.f  () 

if  i:<  .Lt.  . •  >  < u: i 

1£L(2  ,  :  )  =  r’L:‘;  Tc<  u; 

ELSE  ' 

3£L(£, litres--*0-  HE 
E:i3Ic 
C<=F lhF(l 

: f  <;<  .  j.  . 1 1  rnt  j 

3-L(  ?,i)  =ri*cT  r‘r  :• 

else 

0SLf3»i>=^>-=T=  re 

;NDIC 

CK=S :nc  o 

IF  cr  .LE.  )  I UE  4 

DELC-  ,:)  =  re>‘-Tc-  ">ge 

ELSE 

3EL(  •  ,  :)  =  r',3--'-,r.  r.  c 
EnOJF 

-  <=F- ;  *jf  ( ) 

ic  c<  .  i  r»i:  i 

DELI:  ,.)=c**  =  7-  :J 

else 

3  £L  <  »  1)  *P“-=TS'  r  -4 

EN0;r 

C<  =  F  tt.F  (  ) 

if  i:k  .le.  .• )  r «e  * 

OEL<E  ,  ;)  sF  <3*tT=*  ’■■HE 
ELSE 

OtL(if:)=r-,E-r'r'  3US 

E  NnI  F 

’PIN*  ,•  ENI'IAL  -:S-:"A’E  <?f  THE.  ►  EL  ATI V  £  £LE"t.;.T  7  E  3  f  3R  • 

« ,*  v.t'  -.r*T:  as:  i-.’ed  -  i»  -•  r-:jOTc»  • 

PRINT  •  I  •  EEL  *  ff'*  .  (i  ),  "rL  <E  ,1)  ,0  iL(  :  ,1)  ,HEL  (•  »  .  )  f  3£.  (5,  1 1 , 


3? 


con .  in 
Coum  =  C C"J»  :*i 

•  -  A3A..  '■•-•■SJ-  r**-:xT 

?uhi  nzz,  *i»gl,  *-  „‘l,::i£,:o  * 

=>ri>:t  •  ,*  10=  if  .£  l r  it  or  I'iiu.l  s:r  ' 

•  ♦»»■** **■* -**■• -  -  noise  •  cFut-t  . . .  *•*'** 

IF  (NOFfc  .15.  I>  IH£M 
NOSR=  (KANF  ()  -  .  .  >  2‘ 1.7  32  V 

N0S7  *=  tRIF  c(>  -  .i  )  •  2»A.  ’7  V  7  .l  E 

CALL  AQC  *<  ( NCSr.)  '(OF  !\  , »  -NGE  , A 1  £  ) 

ENDIF 

op.IHl  INITIAL  ME*SUti"-Nt  *  ,Ct’l!‘T  ,  *  OA  i  A  ' 

3*I‘.T  •  *'  R»"lGr  s  SiANO',*  c  rl'.Tt.  =  *,Kt.ATC*'  TiM*  1 »  f  I N  F 

=>R1K  r  *  ,  *  • 

^<N=RAf.G£*  >3  ~  1  ,i,j 
^PKMsRLATE*  7.E  I  3;  i  3  '2 
ri«r“  =  '.  I“E*1  •  f-i  :  5 

»M»Jl  ’4MG£<K*>=  »,3  .  Ff,  *  r  I  K£  <  1 1  *.) 

*:  ',7jrE*1 
sto--:=i top  - *i 
*NG(STO~)  =»*.  G£ 

(?KT  (FTCKS)  =  \?A  Tc 

tsto;  (F tor  ;>  =-  im- 

ccp-F"f-:  i n r e - c r- T'* '  -■’ciiiicN  4.0  /elccitv  *• 

L  =L! 

CALL  <t  rLf  (TIUE,LCC,F) 

L S 1=  ■  » L  CC *  ■•  t  ►'  ( ;  ) 

call  Ft  s'/fi. <l  :  ,ci ,  :-cc, 
call  “ixtj  ino',;>  =  o$,  jof:,*> 

CALL  "tXT!  (  '  t  T  I ,  </  £  L  ,  7 r  LT 
00  Z:  1=1,  - 

°ciS<r,r'3-r  ;  =  <: 

Vri  S  II,  5  ID1  •  >*V;  1 1  (I  ,i> 

31  COhTINI". 

1*3  C 0’«T  I  Ml ’  £ 

•*-  co‘Fir-;  *  a  ;  get  3c‘:*'.v'.  -tc  v-iloc:’''  *  • 
if  ( t 0  ,io.  - )  i •* : ►. 

STCKi*ITCF  .  *  * 

ENOIF 

ci.=:i.  *i 

xx=y*  »i 

L*Ll*f*LL(i*t> 

S  =  GI  OIL  13, 1) 

cs*'sio-:l  li, :  ) 

*  =  M10t  l  If  ,  1  ) 

HS=UF  I  ‘  Oi.  L  (3,1) 

£C073T=S0fT  <  »G 
NEAr,‘0  =  (1,/.*  :> 

lst**’*  an-.c-*  rr:  1.  ■  <  *  to  :  >*o;'l  <:,  d 

CALL  FCT:3(=XTT) 

£C£f.T  =  l  37 

50  lscs£C£M-:c::h' 

0LF7‘  £  =  i •::?*  f.c  :  os  i .  .T' 

EC £1  T= LCEK  7*  (L  ST-.CC  )  '  rL  ZC'i 
C<=l£T-L5C 
C  <  =  »  '"S  (Cv- ) 

IF  <C<  .Li,  1 r -1 E  »  G1  *0  f 
30  TO  : 

60  CONTINUE 
E=£C£NT 
LS=LFT 

CALL  FKHT7S.a 

CALL  °C'i/tL(L  *0,  r  CCT  ',7  ,£,*Ci:  ,  JiU  ) 

_  .  CA.U.  M:l,T7  i.Arj7,£;fT  v3-.c_-:.  ■, 
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*•  .  _  i 


■r  vilocity 


call  -L  l;  z  < :.u'- ' , . i ' , ' 

**  c ooc,,T r  *;  r:c:T.v  \ 

oo  '  *=: *  * 

ST  A  7  £  ( I  ,  1  >  =  f- 1_ 3 i  l  (I  ,  I  >  ~  "r  L  C  ‘  C  '  -> 

S"A7-.{i*?,l>=9:»7I(l  »i  )  -1' _ 5  ('  l  i  '■  J'  £  > 

70  CONTINI'i 

»**  KOS  lINcAF.  .SAIT  •-T.-T-TlOr. 

CALL  H!"  <snrs,  Hk‘  1 f  ) 

CALL  MILT  i  C  XX'*  Tf7  >/•!•  Ji  UCL'L  »  ,1  j  -  > 

RAN3-'G=S0PT(S7A7:  ( 1, ->  "  ■  T .  (  _  ,  . )  •  L*STA  i  C  ,  1 )  *  *  2 1 

rra:  -  g=(l>t  ust  ,1  >  s*:  Tr  c , ;)  mc-tc  n  *stat  c  n  .1) » 

♦STAT  £  <  t ,  1 )  -  ST  t  T  £  ( i  ,  i ) )  /-  .,r  f.  £ ' 

00  si  :*i, 

HH (  X X  1 1  )=HC.'-L  (1,1) 

81  OO'PrNU'i 

00  C.Z  l-lt  3 

hm(  xx n ,  ;i  =ijcu;  .  (2,: ) 

82  COHTINI'E 

«♦*  r  C.-1  F  'JT  £  -  f  f.  3 1  1X0  G:  -ATS  -C310UALS  v  "  ' 
R£S(xX,i)  =  •  X3(  L‘,  11  -T'-I’S- 

R£S(XX+i,l)=  T  (i  r  C'-.i)  •  i  .if  .7 

I p  CO  .  in .  1  )  G3  '0  : 

IF  (Cl  •t''.  '.uu';l)  30  "C  !.; 

30  TO  ' 

88  CO!IT!Nt'i 

i  o=; 

RMS-  = 

RMS-- 1- 

oo  :=i 

“S“  '-3(1 , 1)  -2 

RMS  FJ  =  vhv*r.::  (  -  ♦ 

9a  co'Cin"-: 

RH$-  r  SO  <■  T  (  X  V  /£  J  »«C  1) 

RMjr-iso-K-'-v-r/’.'J'' '£ '! 

POlfiT  1  )  *  -  *S  ;  A  %  3  r  CX  2  -X  •;  G1  -Ail  Rl31“i<JAL  £  rv-  L  ■  S  TiJ)R£ 
»KI"  -  13-  »  SU-IS.C,  •  i,s 

cycl  ;scv:i 

•»-  •  C“  n:._  ••  ll«  ;*  nc.f£0  rsc  1- not.  ‘cj-Cls  *“ 

CALL  TA  =CS  (XW,  WH-, c  ,  .<■) 

CALL  K'L'l  >(“•'■.  F,-.V.  I,  "H -c* ,  -1,1  ) 

CALL  klLT]  1  JTr  •  t  *  * 

CALL  LI  x/rc(  M*  -r~,  -  ,  •,  u*  '•'“••a  ,  -  »  K"'.  '  i-  *  1>  ) 

call  pult:  «(mt=''i- <T 5u*p »  »l) 

CALL  rjLT  1  <(  =  •' '  s  t  -TN  -u  'c^  ,  *  -  >  i  ,  3  1>  » **  '  > 

CALC  ‘  M'LTl  »('JC  =%  :  E  'r'-r'L  , -L  ,  1,  > 

PRINT  ,  *  "j:-  .h:  sC-E  J  '0  *u_  3' t  ;  VICTOR  C^XRL't_D  5y  iHi 

call  xCk :  ( ■  ( :  :  l  -  i  ,i 


FI.T-.R 


r  I.  T '  R 


»•»  LiROfTS  1-.L:rTv:  O'TCTtL  :Lt“:  I!  3TATC  VLLi'.'f 


oo  °i  :  =1 ,  - 

OIL  C  ,1)  =T£L  C  ,1)  ♦  -CL"  -.L'  1,1) 

r  L:-L(:,i)  ---r  ( :i.‘  ::.d  > 

91  cominls 

oo  c  j  =1  ,  • 

00  9 c  J^'.i  5 

=  H  (1,  J'  *Hi  C"'-  C  ,  )) 

92  CONTlNi's 

93  CONTI  MU  i 

UFC<T;  3  at-  3  £T 
L 1  Li  * 0 1 L  ( 1 , 1 ) 

G*OI*OtL ( 2 , t ) 

■3S=CS  I  *0£l  (•>,!> 

Hix;.ra  (C ,  i ) 

hs=ut  :  ‘C-.i  r- , . ) 


HEAMC  =  EOF  7(  1  .  /(  L  '  •  ;  )  I 
CALL  rC7o •«?:“) 

XX=-1 

DO  S?  1=1,  :.J ‘■•'EA 
XX=XX*2 

•**  F.iCC*3UT£  TAEGET  FCSTTIGN  ,-t.C  VELOCITY  E.iiC  J?0(  *•* 

•*-  UPDATED  lLEVENVS  -*► 

LST=V£f  MCxTSTOV(  [  )  *o;l(  I,'.) 

EC£N'=LSr 

lsc=£C ent gt»  e:*i  izz~ 

0Lscre=i.--::cTG7  :us<ec- ■  -i 

SCENT* ECEKT*  (L  GT-.SC)  /  CL  IC'i. 

CK=L5T-LEr 
3K=«',S  (CO 

IF  tZK  .IT.  K  -•  >  G  j  TO  9  - 
SO  TO  c... 

COM  !  NU c 

e=ec:m 

call  pcgve l  <l  ,  g,  £  r c*  ,r ,  =os * ,  /«  lt  ) 

CALL  fCLTl  <  7077,  =  E  07,  =  "ST:,  3) 

CALL  MUL'  I  t  ?C  ~  » VI  l  T  ,  V  -  LT.J) 

**•  P;CO“3,JTE  7  EL  A  I  7E  PSiriCH  MO  VELOCITY  STATE  VESTS  fi  ■»* 
00  "7  J  =  1  ,  7 

S7:it  (j,  ii  =  cos::  <  j,i  <  j,  i) 

STATE  <J*  1)  ='/'  .’I  (J  <  J,1  ) 

corn :  vue 

?AriC'EG  =  SCF*(£7;T;  (i,  L)  •*  ?*5T  -  T  c  (  2,  ; )  *  2*  STA I  E  (7, 1 )  •  *  2 1 

RKAT  EG=  <?T  ;T  i  (  , : >  '  S  71  7r  ■! ,  i )  *3 7  j*  i  (.  ,  1 1  «  JTl  T  E  1 2  ,  i )  ♦ 

‘ST  A7;  <‘  ,D  *S* --7( -,D)  7 

EC-FJTE  -.A  VS;  VO  'A"G*  -ME  OSOUALi 
testtx  ,D  s  .vgt;>*o(is:g 
3ES<yX*i|l)=r-7<;>  -•.’’A  Tr  7 
00!.'  HUE 

VHS.  -  =  ■ 

no  07  1=1, 72,; 

r‘:F=  •:=(:,:)•  •  ie~ 

7  r?r  =■  E  ,< I  *  1, :  I  r3 

COT  .  Ml; 

=>M0F  =  S'*  AT  (  -yS'  /r.U*“S  D 

AMC-  -  =T  0  T  (  •  “!  .r  /  'VI-  1;  £) 

’RIVT  *:  -'.a VS;  '!  { ' :  G  E  -,A  T  ;  TE5I0UAL  £Fi  2i.$  *?!;*  FILTER  • 

A':!'.*  ,*  =  Nr  !3‘,  *  EV'f.  =  ',T‘';rrT 

0  0  ‘  ?  E  =  1 ,  '  1 
“.=  o(.  ,i>  =e?:  ('  7;  ii 
C  0N7  E  N 1 1 E 

;ci:l  =  10vT  MT--  <*«!  <;,.)) 
c c  ( 2 )  =  l  c-. ■  (-ir-  "*•:  < o,  :i  i 
CC  (Tl  =  S C  •  T  (AT  F  0 Hi  (3,  -)  1 
CC  ('■  I  =  CO*  T  ( -*7  =  rwI  (  ,  )  ) 

CG  (  I  =  Lr  -  7  (  A  -  (  ■  ,  II 

CC  <•)  =COO  ( H7  r  I  (.,  I 
•  »*  E  L '  V  JG’  •  C:  C'T-  Cy  •  ‘ 

COll!  T  = 

00  ;  1  =  1,=, 

IF  O  ELCEL  (I,  I  I  0r(.l)  THE! 

i out.:  =ccn  7*i 

;’!D1F 
CONTI  HUE 

s^iv*  -  coo  fc<  celt-  check  • 

CnT= 

oo  :?•  i=it32,: 

IF  (,ES(I,l)  .L7.  ?-'IF-)  7U1(. 

:,(T=:n7  *i 
EVOIF 


*  ♦  * 


t  '  - 


•  ccfT'  ii-  t  -G-  - lisijualg 

REC<YX,1»  =■  M G  C  :  ••  •_)  *.r* 

^£«0'X*1,1>**--T(5ir»  •)  •{?•;'• 

IP  (  30  .if'.  1)  GO  TO  ;  " 

IF  (Cl  .£,0.  IU“^1)  GO  TO  ILL- 
GO  TC  32u 
36C  CONTINUE 
RMSk=C 
RHSf>R  = 

00  362  1=1,32,2 

R‘,?R=RE«  (t,--)4  *  2*RIS°. 

PMSFT=.RE3(I  +  i,l>  ►2+-N5r^ 

382  CONTINUE 

*HSo=SCFT  (  '.*$=  /NJ-NEM 
RMSFE  =  FOOT  (?''<'•  P/vUM  |£  .) 

PRINT  vis  RANGE  HO  ’ANG;  R <- T E  P.GI0U1L  EULrS  3Er3?£  FILTER  * 

PRINT  RiS:  =  *,R1;T,  •  -R?KT  =  ♦,  <MCrR 

I01=? 

CYCLE  =CYCL I»1 

PAINT  CYCLE  =  *, CYCLE 

»*♦  EAYEC  EETI-ATUN  c"U*TiCN  *-* 

CALL  T f  PC  S  ( HH ,  PH  1  5  ,  B . ,  ,  ,  L  ) 

CALL  FULY1  \(  HMTF,  OIL  V,  ux-r",'  t  ) 

CALL  PULT  I  .(pn-'TFO,  HH.-I-CCL,  =,  ,r.  ) 

CALL  MULT1A(wuT=,-5IN;,  HI,  =  ) 

00  3  j  1  1  =  1, j 
00  3E-.  J=l,( 

PF  <1,  J)  =P' C,  J>*  HT  =0H  c,  j) 

39?  CONTINUE 


391  CONTINUE 

CALL  LI(i'>2-(cF  ,l»:  jP’l  ,  ,  »*  A-  ££  ,  I  Er.) 

CALL  HULT1  >(•<:,  kz  *,M  T,  *1  , f ,  1 ,  <.  ,  “  .  > 

00  -  I=l,i 

OELOiFll  ,1>='0ElMI,1)-'iE-  t:,i) 

403  CONTINUE 

CALL  M’LTl  !  OR ,  D£  L  SI  F ,  H 
no  •  •  1  I=1,E 

TO’  C  ,1)=N;  C,l)  +M  ,  (  I,  -) 

410  CONTINUE 

call  MULTI  (3  0.  ,T;r,0£.  C£  L,  ) 

ORINT  •',«  THE  CHANGES  TC  r,J--  GTlTE  VECTOu  COMPUTED  3Y  ME  Fi.TER  • 
CALL  MPc:TT(C.LD£.,6,1  ,S) 

UFOATE  N  £L  L  1 1 V  £  0'-cITlL  EL-rENT  STAY;.  VECTOR  *  *  * 

0  0  -  Z.  1=1  ,: 

OEL(I  ,1)=C  :L(i  ,1)  YCELOEL  (I,  1) 

OELOEL  (1, 1)  =  if,EC:LO  EL  (1 ,1)  ) 

423  CONI  ISLE 

ICC)  =50RT  (P'i  (1, 1>> 

CC  (2)  =6  OP  T  I”:  (C,  2) ) 

CC(E)=‘C=T  <’ F:  <J,?)> 

CC  (  )  =SC-T  >P0'.  (•,')) 

CC  (r )  =  SPPT  (oo-.  (£,:)) 

CC  C ) =EOA  T ( 3F  i  U ,  :  ) ) 

00  -  5-  1=1,5 

IF  (OELCTl  (1,1)  CC(r>)  IH<.N 

COUNT  =  C CUN  T ♦! 

ENDIF 

430  CONTINUE 

°kINT  COU1  o  Z  ' ,  IS  UN  ' , '  rcr  Oc  LTm  CHECK  ' 

CN7  = 

***  UPDATE  Tlf'-ET  E.£M-.'TS  *  * 

L  =LI  *0 LL  ( 1 , 1 ) 

G  =  GI  *OEL  ( I,  l) 

G  S  =G  S I  ♦  0  E  L  ( ,  ,:) 

H=HI*OEL(r,l) 

_ |  fi  _ _  _ _ _ _ _ 
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>  S  c.  U  J  =  3 


ECCTr-T=so?T<  ■  .-0’  •  2 /.-••?» 
i£A:.-,0=u0f'  7(1. /(.  ' 

**'  f  SCO-’^UV  -  Tl-r,-'  •  «»-r  '<  iL’jCI'r  ' 

»*  UPDATED  i.Li'-N75  “ 

CALL  RCT33C?C'T) 

XX  =  -1 

30  1.3?  I  =  i,XUf-N£i 
XX=XX ♦£ 

LST=MEANUC>TS1 OP ( 1) + DEL ( 2 , 1 ) 

SCENT  =  LST 

lsc=ecent-ecc:  gt'  -in  cecent) 
olsco£=i.-ecctgt-  :c3(ecent> 

ECENTr  ECENT*  <L  ST-lSC  )/  EL  EOT  E 

CK=L$T-LSC 

CK=AES  <CK) 

IF  CK  .LT.  IE-1  1  GJ  TO  E7f 
GO  TC  ■  3- 
continue 
e  =  ec-:nt 

CALL  POSVFL(L,  G,ECCTGr  ,r  ,'"CT  ,  VFlT  ) 

CALL  MULTX(?CTT,FCST,- OS'T.2) 

CALL  MULT  I  (  TOT  T,VELT  ,  X  EL'l,  3) 

•**  PECO 3 '  J 1  I  EELATl  7?  FDS1T1 UN  ►'ID  VELOCITY  •"  l 

00  1 E6  J=i,  7 

STATE  (J,  U=FOST  I  < J,l)-=rDS(  J*  i) 
ST4TE(J+?»l)=VEwTI(j»l)-7E.S(J»l) 

CONTINUE 

RANC-E  G  =  SOFT  (STATE  (  1,  .)  **  ?*>7ATi  ( 2 ,  1)  •  *  2*ST(  T  i  (3, 1 )  ‘ '  2 ) 

RRAT?G=  (iT  ’.T  •  (>  ,1 1  -S  /A  TL'I,  A)  ♦tiAT  E  {:',i>-S7ATF(2,i)* 

‘STATE  (t  ,;)*Si  ATE  <2  ,D)  rr  ,f,£G 

CC-FOTt  .-CAN S3  tXD  ‘A'JGl  \H  E  JESiOUALS 
RES(XX,i)=  'SG(I)-Ti  "GEO 
3  ES(XX*1,  1)  s'-  -  T  (I  )  -r.  -4  "r  *. 

co.riNUE 
RhSF  =  . 

RMSPr= 

DO  ■  :£  1=1, 32, 2 

•-,"5  R  =  F  ES  ([,!>••  ' 

»r£KF*-Esc:*:,i)  * 2* 

C3N7TNUE 

3  X5‘  =$0  =  1  <■•»$-'  /*''J“-ME  .) 

RM.FF -  =  £0=T TS/VLN  ‘ED 

P«4I*.T  <DS  -AnG’-  .X->  ANGE  -aTc  (E3..DUAL  EFK3..S  » ^  T  E  ?  FILTER 

°<INT  -XS'  =  ‘,{15  7,  '  =  *,-2HSl>R 

CONV:ir--NC:  CHECK  * 

DO  ’a  1=1,31 

Ft'  (I  C-i  S  (I  ,1  )) 

CONTINUE 

00  '  1  =  1 , -2,  c 

IF  (•-£!<;, il  .LT,  T-CIGT)  TH-N 
ci-7=C'.r*i 
E':D  IF 
CONTINUE 
DO  ■  :  .  1  =  2,31,2 

IF  ( :  ES  (1 , 1 )  .LT.  7*c: ':■?’)  TuEs 
CNT  =  CNT  ♦! 
sNNIF 
CONTINUE 

3  RIM  C  DU*  ^  =  ',G>l7,'  =cr  :  E  SI  ’“UAL  CHECK  ' 

IF  (CNT  .ED.  El)  THU  I 
GO  TO  it 
ENDI F 

IF  (CYCLE  • ED.  * )  'HEX 

P'INT  FFUGSA“  ;?:r*  If  GirUif.riAL  “ODE  ' 

GO  TC  «?? 


•A--V 


513 


999 


; jU"i = 

ci  = 

<x=-i 

STor  £=. 

90  T3  ;2 
CONUNUS. 

*♦-  output  fim;.  c.euent  ltat. 

print  *,*  rilati  4t  V.  i"  cur  sthTl  \i:z  tck  is  • 

CALI  MFKIf.T(OcL,  =  ,  :  ,  •>> 

print  »,*  cjvf.Ait'JCf  if  » 

CALL  MFPIT  T ( Pp 1, f  ,«,■•.) 

IF  (10  .50.  9«)  :-0  TO  1  >« 

CALL  LiNVcP(°P:,(,£,PJ  •_  t  ,  i  T  r  > 

00  !  1;  1  =  1,5 
0ELMI,1>=9ELC,1> 

CONI INUt 
COUM= 

C  UT= 

37CFi=i 

CYCLE*' 

Cl'--. 

xx=-i 

50  t:  3r  > 

continue 

PRINT  *,♦  RAi’AR  CAT*  AMI  r:UN  CC  raL£T  £  * 

End 


85 


2ET“i 


Subroutines  Used  in  Computer  Programs 


REAL  l  i  L  5  i  3  t  1 , r  i  ■*  3  •  r  i  C»  “-U 
real  ,;s:i 

CO'-rcf;  /.LrH/  L,  l  4<  ws> z,  ccc ,  iu 

X POM »!,:>  =  ( (L<  Z)  /  -I  J*  •('■>'(  2)  -:CC) 

XP^w  (2, 1)  =  a-rz-ui  •  t  J  MO  ) 
xpow<;,d  = 

XPQW(L  « i)  =  (-MtJl'i  IN  (*>/(',*  (1  .  -ECO'  Z021  £)  >> 
XP0W(5 , 1)  =  I  MU)*  G'  23  =  *ECC*C0S(  "1  >) 

xpnw(6ti)= 

END 


ROTATE  COfPJT-S  t  ;  X RO'A'ION  MATRIX  FOR  THE  POSITION  VELOCITY 
VECTOR  FROM  THE  ROW  Tn  TH"  1J<  F F if  S 

SUE" OUT  INF  ROTATE!-) 

REAL  L  ,LS,3,G3,U, -IS,  2,-00,>'U 
INTFSEI-  1,1 
REAL  RC.:f  ,11'  ) 

COMMON  /:LM/  L, IS, G.  US, r,ECC,NH 

R  (  A,  1)  *C3S  MS)  -  CC  S  (  33  )  -SIN  (HE  !  •  EINIOi)  *  h/G 

R  (1,2  )  =  (-!.'.)<  (CCS  (HJ)»S7*I(GS)*ST).  <H3) -COSt  S3)’ -i /G) 

R  ( 1, 1)  *$i)  MS)  •'SORT  C.  .  2/G‘  *  2) 

R  (2,i)=$IN  M?)  '  CC:  (O'  )  ♦O')*  (HE)  -iIK(3S) *H/G 
R  (2,  CM'TI  ;(HS>-‘  In(  IS)  *'0S(HS) 7  COE  ( OS)  *  H/3 
R(2,3)  -  -CCTfPJ  )  •  CMT<  -1  ,  ->.H  >  2/3-  *2) 

R  (3,1)  =$:  N  (GS)  -  sr>  r  o  .•--(*  'J/6'*  2) 

R  (3, 2)  =  CCJ(3E)*S'RTf  1.  2/G**  2) 


■  <3, - 

)  PH/F 

R  (-■.,. 

)  =  P  ( 3  >  3 ) 

R 

)  =P ( 3, 2) 

R  (  3,1 

)=P (3,1) 

<(£,: 

)=•'(?,  7 1 

R  <  : 

> *R ( 2, 2) 

R  (  >, 

)='  (2,1) 

R  : 

)pR(5,3) 

R  (-,': 

)=f  (1,2) 

R  (•*,  • 

)  =  R  <  1 , 1 ) 

00  2 

2*1,  ’ 

00  1 

J*-,! 

.R  ( I  ,  J  3  = 

I  CONTINUE 
o  out  r  hue 
00  ■-  !*•  ,  i 

00  Z  J  =  1 , 3 
r»(I,J>  =  r 

cor.'  nu: 

COM* 1NUE 
END 


o APT  C0M°UT2S  ‘HF-  NcRIVCTIVE  lAT.\i:*3  Zt  THE  rotation 

HAT-  IX  WITH  RCSPFir  n  TI--  CEl-UNAY  ELEIEnM  G,r,H,(s 

SUBROUTINE  P«  -  T( CRS3,  p:  PH,P*PH$) 


THIS  SUERCITIF  2  C0'3J'TS  P-'Ur  "A  7  3  C C ES  OR  =M’<THLS  OF  THE 
ROTATION  FA  TF IX  «(Tw  •’-S'"''!  TO  G,G5,-t,HS. 


T  »  T:  r-t 1  :  ,  ) 

*L  ot  FC- 1 1  *  r  1 1 !  i'  )  i  3 '  3G~,i#0»l:(  )  !r)»  oR-HS(ll,. 

COM”  ON  /'L  «/  L,  L  >,  3,  V),  J,  <1;,  I,  ICC  ,1'J 


oopG 


5 

13 


15 

21 


PF,°G  13  ,1)  =  'IN  (  WS)  ‘  S 
PRPG(1,2)='TI-(HS)‘: 
sspf.n  .  .21  =  tsi  r:  (HS)  • 


:  i?  •  u/c-‘ 


U|i(-  u1-'-'-  1  -  •'  -  ; 

P RPG ( 1 , 2 )  =  < SI f  ( HS )  ‘  (H  /SI' 

°kPG(2,l)  =  -30S(H?  >  ‘S'  ''(OS)  'H/G  ‘*  2 

PFPG(2,2)=-00S  (HS> '  ’  H/G  2 

Pf<PG<  2, 3)  =  (-COS  (  Pi  >  *f  4  •  *2/3*  •  3)  )  /S'RT  (1. 3-1*  *2/3 “*2) 
PRPG(  3,1)  =  (51  >  (GE  )  -  (1  *’  2 '3  ”))  /aOF  T  (l.--H“ 

prp&(3,2 >  =  cr;  (G; ) ' (h  *?'3'» :> )/sc:r (.., -h*» 2/g**2> 

°RPG(  3 , 3)  =-rt/G-<  2 
03  1  1=1,3 

00  :  J  =  ‘  ,  £ 
pr?G(i,  J)  = ' 
co.t:  nu£ 

CONTINUE 
00  2 

00  1-  J=  1 , 3 

P'f  c  G  ( I ,  J )  = 

CONTINUE 


1(1.1 


-H'*2/G**2) 


CONTI  NOE 

PRPG(T  ,1)  =  °RPC-  (l,lt 
o  oPG  ( f  )  =32°G  (2 , 1  > 
PRPGtE ,T )=PR=G(3,: ) 
PpPGd  ,‘  )  =JTPG  (1,2) 
=  F  OGC-  , C >  =PRPG  <2,  I) 
PKPG(b,N)=oRPG(3,2) 
PRPGtl,,>)  =  3’cG(l,?) 

PRPG(v,t)='TPG(?,3) 
Pl.OG  (E  ,(  )  =  >RoG  <2,  ?) 


OP,  PH 

OPPH(i,l)r(-£II(H5)';','<("?'/G) 
prph(i,2)  =  (-r:\(Hj)  ‘:o-('-5i /:•) 
a  ROW  ( 1 ,2)=(-riN(H5)*-,/T»?)/SC'.T(l..-H'“2/G“2> 
PRPH(2,1)=C0S  (HI)  '?*(  'IT' 

PRPM  (2 , 2)  =  C3 S  (  HJ  >  *  C0:  '331 

opcH(2,3)  =  <cc:-  (HP)  •  .■>)  /cr  Jf  (1,  .w»2/G*  ‘  2) 

0RPH(3,i)  =  (*SI  )( C  5)  *  J  '3*  ’  ?'/ECTT  ( 1  •  -  -H'"*  2/3*  ‘2) 
PRPH(3,2)  =  (-Coa(3-)‘T/3*  ?>  /2  0RT  ( 1 .  J -H»*  2/3*  *  2> 
_  o ROM (2, 2) *1.  /G 
00  2;-  1  =  1,3 
0  0  21  J= 

P'PH(I,J)  = 

21  COtriNUi 
25  CONTINUE 

00  T:  1  =  -  ,  O 

00  3  J=l,3 

PP0H<I,J)=' 

31  CONTINUE 
35  CONTINUE 

PRPH(w,i  )=PRPU (1,1 ) 

PROH(E  ,>  )  =°RpH  (2,  1  ) 

°R°H(  £  ,  *. )  =P?PH  (  3  ,  1 ) 

PkPH(t,t)=P?PH(l,2) 

PRPH(=  ,*)  s=RP»  (2,2) 

P).PH<l,‘.)=oRF*MZ,2> 

PRPH(L  ,>,)  iORPH  (1  ,  •  ) 

OPPH(2  ,1  )  ='JTPHC,  ;  ) 

3PPH(E  ,b)  =p?  =  u  (:,  :> 


,1« 
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?rv= h: 


!  4" 
|45 


St, 

55 


63 

65 


7i 

75 


3  :  pM  3  (  1  , 

;) 

=-  ( .  :*•  ( -< 

°  RPM;  (-.  , 

C  ) 

=  i!  ■  (MS  ) 

P  P.PH  ;  ( j.  f 

7  ) 

=  ;■ :  (mm 

Of.  PW s  (  2  , 

1  ) 

=  Z0;  (Hi.  } 

PC;PMS  (2, 

2  ) 

=  -  (  C  C.C  (  M 

PRPMS (C , 

3 ) 

=31 ‘ (MS) 

—  )*<'-/3>'S:M(GS)‘3  0SM'')> 
)-C23CK>  J05(GS>--l/O 
.  -KM2/G-*?) 

>  -i".  M(  MS)  •  31  I  (OS)  *  -I  /G 
'-M  *S.M  («i)  *CS»  (G5>  ‘H/-,  ) 

.  -u-  2/G  *3) 


00  i  y  1  =  1,3 


00  4.  J = 4  , o 
PRFHSCI,  J)  = 
CONTINUE 


CONTINUE 


00  13  :=(:,= 

00  5  J=1  ,  3 

PR  =  HS  (1  ,  J)  = 

continue 


CONTINUE 


®rphs  (3,1 )  =: 

:>K=H£  13,2)  =  : 

PF=M5 (3,31 =. 

=>RP“3  (<  ,s  1  =?  =  '  m£  (i,H 
PPPM3  (:  ,.  )  =Pr.  =  HS  (1,11 
PRFMJ  (t  >=°FrHS  (1,  H 
PRFM;  (4  ,.  )  ==>-PwG  (i  ,  ?l 
PF-PP!  (f  iPppuj  ( ■>  ,  ?i 

PiPH3(t,f)  =  P=FHS(3,3' 
PRPHi  (.  ,G)  =  Ir  =u£  (1 .  7> 
Pi-Phi  (:  1  =or  =  MS  (i,  ’I 

PRPHi  tt  »5)  sPPPHJ  (*,  •»! 


=>RPG3 


=  I  CG3  (1,1)  =  -  (003  ( -t>)  - 't '((*-')♦$;(  (h£)  -  CCS  (G3>  *m/-.  > 
°P?  G?  (.  ,2)  =31  •  (MS)  -  O'  ‘l  (03)  ■  H/G-C  0$  (MS  >*  COSOS) 
pP.=GS  (1,3)  =  ' 

pP=C-S  <  2,1)  =-SI  MM3)  43  *'(  (33)  +3  OS  (MS)  ‘COSt  GS>  ’  H/3 
pp.PGS  (2,2)=-  C  :i.  (-*3)  *  ,'03i  3h)  +C0£  (MS)  -  SIN  <G3'  ‘M/-,) 
PM.PG5  (2,3)  =. 

3 PPG?  (2  ,  i)=CCC  (Gi)  *  pi -’M  1 .  '  2/G* • ?> 

pp.PGS  (3,?)  =  -$.MGi  )  ‘5  0,T(  ‘  ,1  -h-»2/C-*-2) 

PRFC-i  (3,3)  =  ' 

00  65  1=1, 3 
00  (  J=  :  ,  3 
P-'PGS  (I,  J)  = 

0 0(27 ;  M's 

00VTI.U2 
00  7;.  I=-  ,‘j 

00  "  J  =  1  ,  7 
p;  PGS  (I  ,  )lr 
cont  :  muc 
CCNTINUt 

P p CC- i  (•  ,.  )  =p=  p  C-S  (1,1) 

P  P  PC-  3  (  S  ,(  )  =  P  7  p  GS  (  2 , 1 1 
PR=G3  (».  ,  >  )  =PP  ?  C-S  (?,'.i 
0  R  c  0, 3  ( F  ,3)  =  PRr  CS  ( 1 ,  =1 
=  F,PGS  (t  ,  )  =Pi  FC-S  (I ,  ?i 
°P°C,3  (t  ,r)=a  =  PCS  (:,7l 
PP.PC-3  (1  ,t  )  =Pp  =CS  (1,31 
pRPCS  (5  ,6>s»p  =  GS  (2,  p) 

PPPC-S  (*-  ,f  )  =pp c  C-S  ( 3 ,  ’I 
Sod 


1 5.S  : of_pj.it z 5  m : 2 '  T  *  >( :  >:  position .  a*o  -islooity  vs 
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AD-A1I1  107  AIR  FORCE  INST  OF  TECH  RRIOHT-PATTERSON  AFB  OH  SCHOO— ETC  F/B  2B/S  , 
DESION  Of  AN  ORBITAL  ELEMENT  ESTIMATOR  USINB  RELATIVE  MOTION  OA— CTC(U) 
DEC  01  J  F  ANTHONY 

UNCLASSIFIED  Af IT/§A/AA/ilD-l  NL 


with  respect  "o  t -i £  -r-r'icvi-ouiv 

SU3,  GU'iTM  *r..  (CX5:E.  LXUfit  , C/POE, OVO'lr,r'V.i''E) 

PEAl  ECC,H’J,L  ,  L,C,  0X~  V.'TpCE  (iyiT.,  D  /POi,0 ' DCE,  '"/WOE  ,.3 
COMNON  /Ll£«/  L,LS,  G,G?t  -1 , -‘S  ,  £ ,  tCC  , -M 
DXPP: =-  (L*  ■*?/•  U) '  iMf  r> 

OXOCE  =  L'  G*  COS  (  E)  /WJ 
0  XWC  -.= » 

OVPOE=  (KU*  (SH  (El**2«  *!.*E0"-MU‘  COb(c)  *L»  (1.  '-EC:*C03(i> ) )  / 
*(  a'  (i..  -El-.C"  COS  (:)>'  '  *?» 

OVQOE=-  f  ( L-»  L *  <1.  -si:  “GD*  (")  )  )*  (i'.U‘  5*3IN(t)>* 

MMU*  3*C0S(s>  •  SIN(E>  *_ )/<  <L*  L*  (I.*;  -ESC ‘COS  '  E)  l  >  “  2) 
OVWCEM 
END 


XOL  SOHFUT  E;  THE  3EV7YT TVS  OF  1  HE  POSITION  AND  NsLDJtTf  VECTOR 
WITH  RESPECT  TO  THE  OElAIUSY  ELc-UNT  L 

SUDS  OUT  IN  E  XOL  <OX  =  DE.  Wi  ?  ,  CXWC  "  ,  Dtf°D  E,  DVGD  E  *  PW  '’E ,  3X0. ) 

REAL  DXPOE  » DXPCE  »  DXW3  ",  3  /o^t » lOVOCE ,  O'/WO;  ,  ESS  ,  HJ,  tp?.,  X3P_  , 

‘XWPL,  VFRL,  VOFL  ,  VW^.  ,  0  "It  ,  L  ,  US 
REAL  DXDL <l«s,iSll 

CGNWCN  /SLEW/  L,l>,3.  G3, W,HS,£,ECO,H'.' 

DEDL  =  (<1/EC0)-G*3'EI'J  (E)  ‘L'*3)/(l.  -ECO*  COE  f  E)> 

X PPL  -  2 1  L‘ CCS (E )/l J-r  L* etc/mu-l* L/ECC/MU-G-S/L** ’ 

XOPL  =  G*SIN(£ )  /-111 

XWPL=.  r  • 

VPPL  =  h|j‘  S I  1(E)  C.  ;-rT0‘T'>P  <i)-L/ECD*G*G*CDS  (E>  'L“  3)/ 

*(L-  <1,-  -sconces  (c  I )  )  •  ‘2 

VQCL  =  -l  U'  G  '■  0  0  S  ( E ) '  <  2 '  L*  (1  .  ' -ECC'CCS  (E) ) -L*.  /  ECO* G*C*  0  OS  (E)/ 
*L**3)  /  <L  L-  (l.  -E::«'TKr"  )•*■£ 

1/  WDL  = 

oxol (i.doxpm- r-D.^xas. 

OXOL  (  2, 1 )  =  iXOPE--  C E0L‘  X'13. 

OXCL(3,i)= 

OX'llli.ilPjV':;'  CE0L‘''->°'_ 

QXDL  (£,!)=■  TV  CO  E’PEDL‘ Hi 
DXDl  ( t ,  1)  =  • 

ENO 


XOLS  COMPUTES  THE  0E5TW»VR  OF  THE.  POSITION  AND  VE.3CITY  VECTOR 
WITH  RESPECT  TC  r  TETE-'JNAY  EL  EH EM  LS 

SUER  CUT  INF  XOL  S(  CXP<3:  ,m  : ,  C  XWCs  ,  CVPOs,  DVD  0  E ,  DJ  ''0; ,  D<  DL  S ) 

REAL  CXPOE,OX"DE,  3<V>  ",  3  r^t ,  DVCDE,  0  WOE ,  EOT  ,  MJ,  ‘‘ED-S ,  . ,  3 ,  E,  .  S 
REAL  DXDLS  (1U  ,1S1) 

COMi' ON  /ELS’'/  L,LS,G, GG,W,WS,E,ECC,NU 
OEDLS«i/(l.S -ECC-  COS<  ED 
0X013 (i,l>=OXFDS*OED.P 
DX0LS(2,l)=r’XOC£‘3El.T 
o xcl 3(3,:)-: 

OXPLS  (h  jDpDVpOE’DED.  G 
OXDLSU  ,1)=0VPCE“0ED.'' 
oxoL3(e,i)  =  :. 

END 


XOG  CONFUTES  THE  :LU"U'r  OF  THE  POSITION  AliO  VFLOCirr  VECTOR 
WITH  PEScEOT  TO  i WE  y-LAIHTY  ELEMENT  G 

SUV  OUTIk’E  XCG  (CYPOE,  pxpi r ,  C-x WO E ,r v»0 E , D  V OO E  ,  CVW  nE,OXOS) 

?  E.AJ--  *  c-cr- }  IV-j Jl  y  F  s  S  ,  /  P  ’  G  .  - .  1 ,  E  C  ;  ,  E  ,  L  ,  G  ,  L  S  ,  C E  0  G  ,  ~  <  PO  E ,  0  X  0  3  E  _ 
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OX'.ilEjCVp'oe,?  VCCi, 

REAL  0X03  (ISr  .  HI) 

30“', ON  /-:t  -  j/  l, Li, 3,  ,£,ECC,HM 

osog=  (l ./ego* o  s i v< s>  ml.  l>  ) /(;  ,-elc-cos  t e»  i 

xp<=r=G/ (‘ u  ?c<  )  _ 

X0PG=L‘  SOME)  /-'U 

VP°G  =  <  hU*  STMf  i)  '  C  03<  M  *3/:  /  £ro  /((.*<  1  .*  >£3  C  *05  (z)  *  )  *‘2> 

VO=G=  (L'LMl  .  -SOC -C)S(E>  W;U*CGS(E)  -HU' COS <  £)'  3  (  E  J  *3*3/ ECO)  / 

♦  (  (L‘  L-  (1,1-tCC  •  CC 5  € ~>  )  >  *  ’> 

oxog(i,d=oxpc-*  c:03‘v',3G 

0X00(2,1)  OXOCE*  C103 *■  TO3'- 
0  XOG (2,1)  = 

DXOG(L  ,1)*')\»P'H‘  DEO>V->o-, 

oxogis ,1) =ovnc:-  os:  o*  vo3'. 

QXOG(E ,1)= 

SNO 


IOL  CCPPU7SS  THE  OE?TVt'rVE  CP  THE  INERTIAL  °05TTI0'(  1N0  VELOCITY 

vector  w;th  reef?!!  *i  ru*r  ci  l  a un». y  HLt-.tuT  l 

SUB-CUl ONE  IOL (p ,0<0. ,YT3L) 

REAL  R(lt£,Ul),CXOL(HS,  1*1),XIDL(1«S,1*1) 

CALL  HULTI  (C,CXOL,  <P'.,  S> 

£1(0 


IDLE  CLNPL'Vi 3  PH?  0"’*/X*T'/£  OF  T Hi  INERT;!'.  POSITION  \ >1 3  VE.OOITY 
V  EC1  0 F  MI T H  RESPECT  p  0  THE  CSLAUHAY  ELE-ECT  LS 

SUE' OUTIfl?  irLS(i,0OL0»'<’'">LE> 

REAL  F  (i  U ,  1 1  ),LOL'  («.*••,  S(i)  , XICLS  (n6,H<  ) 

CALL  HULTI(;.,CXCl;,Xf  OlS,S) 
i  NO 


ICG  C  0V3U  7ES  7H:  Di*  ’VI*  T  »£  OF  THE  INERT  It L  poSTtion  tNO  VELOCITY 

vector.  with  tiiFLor  ro  r-t-  celaumay  ELEHcNr  g 

SUP!  0U7INE  ICG  (FF  =  0 ,  <  33W,  5  ,  CXOC- ,  XIU3 ) 

I  NT  I G  £  K  I  . 

REAL  PKFG(1«G,1*S)  ,?(  Hi,'  *6)  ,XFCW(i  tt,lSl) ,  0X03  (Hi,  HI)  , 

*  X I C  G  ( 1 » 6 , 1  s  1 )  ,  1 1 : 1  S  p  ( • ' ) £ ,  *. !  I ) 

CALL  ML'Li  1  (PR3G,X--'W,  T'lr_  =  ,6) 

CALL  rULTI(R,LXCO,<nG,F) 

o  c  i  ;  =i , e 

XIOG(I,1)=XIOG(I,11  *INT£^n,i) 

CONTINUE 

SNO 


IOGS  CONFUTES  7HI  OI’TYH'M/E  OF  7H0  INERTIAL  FOS'TION  &N0  VE.OCITY 
VECTOR  WITH  FSSPfcCr  3  DELAUNAY  ELEMENT  GS 

SUBROUTINE  IDGS(FR=3*  , X30 N , XI DGS) 

REAL  Pr  PCS  <1*1  ,1  S?)  ,'T>M(1  !t  ,isi)  ,XICGS(lti.  HI) 

CALL  HULTKPBFGa, <  =  ■)■<, TT7TG,L) 

£NP 


I  OH  CCHPUIES  The  OE:-w:  he  OF  THE  HERTHL  POSITION  tNO  VE.OCITY 
VECTOR  WITH  RESPECT  p0  Thp  DELAUNAY  ELEMENT  H 


S’J5‘  CUT  INF 

PEM.  F.  bm(i  j;  ,:!;i,ri4i:  ; '.  ,1 :  i ) ,  x:  an  ( : « &,«  •.  i) 
Call  mulu  (c  vc h, x  =  h.  vro-i,'  j 
if  0 


IOHS  COMPUTES  THE  3E’;tfT'TV£  OF  THE  INEPT  IV.  F03TTIJN  4N3  VELOCITY 
VEC’CR  kith  \LSci:r  *0  r-t  =  DELAUNAY  ELEMENT  ns 

S  USE  OUT  INF  IG‘n$(F<045  ,  <0-1 «,  XIDHS) 

PEAL  PPPHSCltr  ,HS>  T£.  ,lti)  ,XI0HS(1«S.  Hi) 

CALL  MULTI  (P«  PHS ,  <=>1  ■» ,  *nn>  ,  £  ) 

ENO 


INEtT  POTATES  POSITION  1  \'0  VELOCITY  VECTOR  cn01  POrf  TO  I  JK  FPAME 

SU300UT1NE  I  SEPT  (  P,  <  vl>)) 

PEAL  R(liC,l!i  )  ,Y  =>'><(  11?.  t  ’1)  ,XlN(i»S,l«L) 

CALL  MULTI(P,YFCW,»T-1,M 
ENO 


FILL  FCR4S  THE  STATE  ’">*  'l?T  Ti  Ot«/COCP',1NATE  T  PAN?  rOPMA  T ION  MATRIX 

using  the  ,,  ocPiv:rrfr  vetoes  cf  th.  in ekt f  sl  3  osi t i d  m  and 

VELOCITY  NITH  LES^TT  tj  THI  GFLAUNAY.  ELEMENTS 

S  UPF  CUT  IN  F  FILL(x:3.,  v,T-'‘,x:DG,x:CGS,XI'TH,<:rMP.TPiNS) 

PEAL  XIOL  <1 1' .,1111  'T'.STi  *0,  Hi) ,  XI 05  (it  6,  1  »i),  VIOS  S  ( 1 »  5 , 1 » 1 )  , 

*XIPH{1SI  ,  if  1)  ,  XIC-*Sf  1  *3,  1  *  1)  ,T.  A  NS  (1  Si,  ISO 
IN'iC-EF  : 

00  i  1*1,  ? 

t-  ans  a, i»=x:oi  <i,i  > 
io  Conti  r,uc 
oo  z  :=i , e 

TEANS<I,3>=XIPL>CI, 1) 

20  CONTINUE 

00  ?  1  =  1, 5 

TF.iN'S  (i  ,  T)=XIOC-(I,'  ) 

30  CONTINUE 

00  -  1  =  1,  *5 

_  T~.iN?  (I,4)=x:cc-S(l,  i) 

AO  CONTINUE 

OC  •  1=1,: 

Tf  *  ns  ti,->=x:oH<:  ) 

50  CONTINUE 

00  :  *=!,*: 

TF  A  NS  (I =x:0M3  (I,  11 
*a  CONTINUE 
EflO 


“OSVIL  COMMUTES  A  3OIITTO‘l  AND  A  VELOCITY  VEC’OP  IN  f ME  3 QW  FEAME 

SUIOUTif.E  P0C  VSLC.S’  ECCCAT,EC1T,P03.  7E_> 

PEAL  LSAr  ,  55  AT  ,EC:S1',”’T,‘U 
PEAl  P0F<H3,1  fl),V*.  (1*T,  HI) 

MU=1. 

POS(l,l)=(  C.SAT'  *  E)  fa'ij  *  (“OS(ESfcT)  -ECCS4T) 

POSCO,  1)  =  <L5AT  GSU/-'  ))  *  ('■TN(EtAT) ) 

POSH  »i)  =  s. 

V*LC .T >  =  <-1U)- -FT V(  =  ;«») /() SlT  (l.1  -SCCSAT'TOS(rSlT) ) I 
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tftL(2,:>=  (  'J)  -  GS-T-3T  -(-'«T>  /(Li  A  7  *r*  (i.  .-ETCS'  '*CJS(ES1T  )  )  ) 
ENP 


RC‘r?J  COMFJT -  b  4  •<*  ROTV'-’-Of  Mt TEIX  c OR  KDTATXJM  FOl  TME  POW 
Fc.tr'-:  TO  IMF  :  JK  FRV  F 

SURT.O'ITInE  R0T3J(OTi 
REAL  FCT  (1  S3,  U3) 

REAL  L,LS,G,GE  ,u,-,3,r,:",M'J 

C 0“'>C N  /in'4/  L,L  5,  ">,  GR ,  J,uS»c,rCC,iU 

R(.T(1,  I)  =  f  ?S(MS)*  30SC  '-.SI-FTMHS)-  SH(G5)*H/T 

ROTC  ,2)  =  (-1.  )‘  (;33'mS)-ftn(GS)*SIMH$)'C03  (G£) ‘H/3) 

ROTC  ,  3)  =5  *S  (u'i)-  SIR'  (1.  -u*  i/O 2) 

RO7(2,l)  =  Sl>t(M0)-  :OS(  '-3)  +'•'$(  HE  )*SJ|J(3S)*H/T 
ROT  (7  ,  2)  =-:ll  (  WS)  •  31 1  (*,$)  i'CE  <  Mf )  ■  COS  (GS  >■*■(/ 6 
ROT(2,3)s-:3S(MS)  >  SO’Td.'-H  -  2/G*  2) 

ROT (3,1  JsSIHCJS)*  SORT  (1.  '  -m*  2/0-  2) 

ROT  ( J  »  2) =C  OS ( GS)  •  SIR'  (1.  -m*  •  2/G-  *  2) 

R07(?, *>=H/3 
END 


1H  COMPUTES  'ME  M  1VRI*  c-.CP  :.£LAUVE  POSITION  'MO  VE.OOITY  CA7A 

shopcutine  MM(sr,Misn 

REAL  PANG (0,1.  F4TL0 

REAL  ST  <1  »:,itl)  ,  -MAT  fit’,  Iti) 

RAHC-EC  =  SO'T(CT  (1 , 1 )  •’  R*'r  f  2 , 1 )  *  2*ST ( 3,1 >»' R  ) 

RfiATEC  =  (3T(i,;)'  ST  (i.  4  >  »t'  (2, 1>*  SI  (5,1)  *ST(E  ,1)  ‘  ^T  (S,l) )  /  RAM3EC 
hma;  (i, i>- sc (i,i) ' r ' -r. f o 

HMA" (i,2)=3T (2,1) ARAAGFR 

hma:  (1,3)«ST(3,1)ARHGfo 
H'  6'  (2,- .)SH«AT  (1,1) 

HMAT  (2,E)=mha  *  (i,  5) 

H  MAT  ( 2 , 6)  =  ■•M AT  (1,5) 

H MAT  (i,<0* 

MM  AT  (1,0=  • 

M KA1  ( 1,6)  s  : 

MMAT(2,1)  =  (RA  IGEC'STS  .,1)  -RK.A  TEC  *  ST  (1 ,1)  )✓(  R AN3* '>*- 2) 

HMAT  (2,  2)  =  (RA*'GeC  *  5  T(  %  1>  -’PATEC-ST  (2,1)  )  /(  RANGE  p*»2) 

HMAT(  2,3>  =  (RAf,GEC»5T(  1,1)  -RRAHC*  ST  (3,1)  )  /  (  »  ;i.3T  V*  2) 

ENO 


°HT  COMPUTE  THE  STATE  TRA MSI 'IOT-/C jOROIM AT E  TFAYSFoRMHION  mATEIX 
IT  A  GIVEU  -EASU= EME(TS  'tM£ 

SUBROUTINE  PHI (TRAMS) 

REAL  L,LS,6,GS,M,MS,E,'rO0,  "M, XFFL, XQFL, XWPL,  VC=L  ,  V  0s  „  , )/  WPL  ,  OE  0:. , 

•Teas,r.xpo',croot,ox'CE,">'/pcE,C'vcr.E,3;w'j£,x5=G,y''PG,x=p3, vop  ,  r£DG 

IMTEGEF  I,  J 

REAL  RdSf  ,1»:  )  ,Fi=G(  tSG,  1  to)  ,FFPGS(  1«  J,U>)  , 

»ORPH(i  SF,i  !«  )  ,FFF-*S(!  !5,  1  s  S) ,  XFOW(UE,lSl),OxrL(H6,lSl)  , 

‘Oxol?  (i  tc ,  mi  ,rx;sc  «s , i »*. ) , im eru  !■,,! m .  xiy.  (i«s, u i>  , 

-XI OLE  (1  t-  ,  ;!i)  ,x::GC  SA,  1  »1>,X1CGS(*  :  ,,i  *  1)  ,  XIE-t*  1 1«  5 ,  HI)  , 
*XIOH(lSt  ,  ill)  ,'F.f  bSC  t »  ,  l !  r  ) 

COMM  or  /EL  EM/  l, IS,',,  *,T,M,MS,E,ECC,i'J 

THIS  SUEFO’JTI.iE  COIPJ'^S  'HE  STATE  TRANSITION  JFo.mritfE  MATCIX 

.ciu.fct.'J;LD_ .  . . 
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mum 


mm 


CALL  *,7  (  »CCHC  ) 

CALL  RV=''*»(*  =  “W) 

CALL  XCE(DX=C  T.CXCCE,  ^XM-I-.CVPU  ,C'VTC-;,OtfW-i-) 

call  xt  l<cx?c  ■  ,cx:oe.  .rvio-:, jvot  ,ov>l> 

CALL  XLLS(OX=~-,(  X ix  CCL,GV<n  If  "I»"'LS) 

CALL  XOG<  Cx  =  0'  (rp?-,  T'-'T'- ,  overt  ,CVOC£,TVMCr  .ex'*  C) 

CALL  lCL(f  ,0y:L,Xi0L> 

CALL  IDLSI'fCxn-Li,  XIT_5> 

CALL  ltGir’er,,  XFC,*,R,  oyV.,vlC0> 

CALL  IDGS  (3RcGS,X30N»  v  r  OG  =  1 
CALL  10H(FRPMf  XFC/,  XT  ">-» 

CALL  ILHS<°^PMS,XeCA.  *IH5) 

CALL  FILL  (XI  CL,  Xi  ?L5,  XTT. ,  xiOGS,  XIOH,  XIOHS,  TRA^S  > 

=TflO 

S  IJPR CUT  IN  t  ACOMNCSR,  ‘OF’-  ,R^  AR) 

SUBROUTINE  SCON  aC0S  TO  M  £  ASUR  EHENTS  S> 

REAL  NOSR.'OST  P,F1,R>  R  <>• 

PA=kA+NCS?  r? 

PAf  =  r  AR+NOS  LR  ‘‘  ’ 

END 


H  ATxl  X- VEC7  OK  MJL  C I  3.  T  OS  *  T  ON  SUERCUTIHS 

SUBPCUTIUE  NULTI  (S  ,0.  b’OI,  N) 

INTEGER  I,  J,X,N 

REAL  ACin-,ltt),r(lS  ,  1  x  ?  >  , PROD ( i * N,  i tl> 
J  =  1 

00  2  3=1, N 

PPODd,  J)=' 

OC  1 .  K=l, ( 

Pr  00  (I,  J)=P  00  (I,  Ji  ♦'  (T,x>  •  3(X,  j) 

10  CONTINUE 
ZO  CONTINUE 
EUO 


MATRIX-NfiTEIX  HUL( I». TTirr^N  SUBROUTINE 

SUBROUTINE  HULT1  (l  4  ,  »  0,  »■>  301,  H,  N ,  OH) 

INTEGER  I,  J,  <,",(■,  OH 

REAL  AA  (1(N,1SH)  ,33  (<  1 X  , :  OH)  , Ph  ODl  (  1 !  M ,  1 1 0 X  ) 
00  I  J =1 , OH 
00  f  1  =  1,  N 
PKCOUI,  J)='  . 

00  1  .  K=1,H 

P-0D1  (I,  JMFFOCl  <1,  J)*‘  l  (i,K)  ’  E“  (K,  J) 

10  CONTINUE 
Zil  C  0NT I NU  E 

30  CONTINUE 
END 


SUET  CUTINE  KEFLEF.  usr~  <E  Pl  FR  A.GORITHH  tc  rf\|o  the 

ECCENTf IC  TNOtALY  j=  AN  ’Jill  GIVEN  A  PARTICULAR  THE  ’AST 
PERIGEE  OF  TH=  OF  3 1  T. 

SUBROUTINE  <E  =  LE'  (II4,T";  ,ECEnT) 

REAL  l, IS, G, OS, H, H 3,*, ICO, "U, TIM, -CENT, M, HI.  DuOE  ,TE3T 
COhHGN  /*L  -y/  L,L5,G, -,3,  (,«S,E|ECC,NU 
1=  (XU*'  2/L  *3)  •  TIr 
ECENTrN 

10  Xi  =  ECcM-E.OOl‘  SIMECENM 

Ohdei1,.-cC1’ plF(EC‘V*t 

-‘Ct*  r  =  £1  r> 


t  _  . 

: *  ■-  r  -  u  --  -•  - f 
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T  SS*  =  - 
Tc'CTiAf  EU  ;31) 

I f  ( .  eet  .u. 

50  TO  1 
20  5oi.t:nu- 
ENO 


MP'vINT  F-InTS  L  PAIVX  O'1  '<ECTQF 

SUOPOUlINf  SF'lNl  O.’^FX) 

INTFSEF  FF  ,5,  I|J, 

=EAl  aU'MAX,1) 

00  1-  1  =  1,  \? 

P  Vi  N 1  *  (’K.iEli.i)  *  ,  (A  (*,  J>  ,  J  =  1,C) 
CONTINUE 
ppiur  » ,  •  • 
eF.IsT  J,«  • 

END 


TRPCS  CCT'F'jrES  THE  INT'nOE  OF  A  MATUX 

308'-DUT:!!E  TPPCSMAT.  11  T~  ,  "’0*  ,CCL  , S  "UX  ,  CMAX » 
INTEOEt  RC’<fCOL,I,  J, »  HX,  '«AX 

?EAL  P*T(F.  1AX,1)  ,1ATT 
OO  l  I=t,OW 
00  1  J=J. ,331 

kATT<J,I)=«ATU,J) 

COH'INUE 

ENO 


1  *  . 

\¥i 


mat^ix-^6 vux  puLr:p.io»*roN  subroutine 

S U3F.0UT  INe  HUL'tlAni,  ,M,N,  OH  ,  NMAX  ,“H  AX) 

INTt  3EF  X,  J,X,  *,(■!,  OH.  I-AX.xmAX 

^£AL  At  (NMK»1  J  ,63 (K  \x,<. )  , PE0C1  (KA4X,1) 

00  ?  J=i,OH 
00  2  I*i»H 
P'OOi  (I,  J)  =  . 

30  1  K=1,V 

PPOOl  CX,  J)=PPOOHI.  M  *’8  (I,K)*BE(<,  J) 

CONTINUE 

CONTINUE 

CONTINUE 

END 


SUP‘  ;UT  IN?  NUITXA  t  A,1 N..#x> 
INTi-GtF  1,  J,K,f„F  «AX 

REAL  A(NVtX,l)  ,  =”00  t  NF  AX ,  1 ) 

J  =  1 

oo  2  :=i,'i 
ppco  (i ,  jy=’. 

00  1  K=l,  I 

P£  00(1,  J)  =CF.OD(t,  J*  »!(*,<)  *C(K,J) 
CONTINUE 


VITA 
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flight  test  engineer  on  NC-141,  NC-130,  RF4-C,  and  UH-1H 
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